########################################
# SPI 
########################################
data=read.csv(file.choose())
#fix(data)
pre1=data[,-1]
pre=as.matrix(pre1)
pre=as.numeric(pre)
pre
#Make 3 month time scale series
library(zoo)
sum3=rollapply(pre, 3, sum) # 3 time month scale 
sum6=rollapply(pre, 6, sum)
sum9=rollapply(pre, 9, sum)
sum12=rollapply(pre, 12, sum)
sum24=rollapply(pre, 24, sum)
sum48=rollapply(pre, 48, sum)


O2=rep(c(0),each=2)
O5=rep(c(0),each=5)
O8=rep(c(0),each=8)
O11=rep(c(0),each=11)
O23=rep(c(0),each=23)
O47=rep(c(0),each=47)

pre3=c(O2,sum3)
pre6=c(O5,sum6)
pre9=c(O8,sum9)
pre12=c(O11,sum12)
pre24=c(O23,sum24)
pre48=c(O47,sum48)

DFO=data.frame(pre,pre3,pre6,pre9,pre12,pre24,pre48)
DFO

length(pre)
#SPI3
library("fitdistrplus")
library("propagate")
FF1=fitDistr(pre3[3:348], nbin = 100, weights = FALSE,verbose = TRUE,  brute = c("fast", "slow"), 
            plot = c("hist", "qq"),  distsel = NULL)

FF1 #GND #bic -829
Log pearson(3P) is appropopriate #bic 3896 #aic 3885

library(stats)
#Fit1=dlnorm(pre3, meanlog =4.1341, sdlog =1.0654, log = FALSE)
#CFit1=plnorm(pre3, meanlog =4.1341, sdlog =1.0654, lower.tail = TRUE, log.p = FALSE)
plot(pre3,CFit1,xlab="values",ylab="Commulative Probablity")

# 2P Exponetial showing missing values
library(tolerance)
#Fit1=d2exp(pre3, rate = 1/0.01838, shift =0.61, log = FALSE)
#CFit1=p2exp(pre3, rate = 1/0.01838, shift =0.61, lower.tail = TRUE, log.p = FALSE)

#4p Burrr distribution
library(actuar)

#Fit1=dburr(pre3, shape1=994.69, shape2=1.0053, rate =0.60942, scale =52811.0,
      log = FALSE)
#CFit1=pburr(pre3, shape1=994.69, shape2=1.0053, rate =0.60942, scale =52811.0,
      lower.tail = TRUE, log.p = FALSE)
plot(pre3,CFit1,xlab="values",ylab="Commulative Probablity")

#Fatigue Life 3p
library(extraDistr)
Fit1=dfatigue(pre3, alpha=0.91399, beta =82.187, mu =-4.0011, log = FALSE)
CFit1=pfatigue(pre3, alpha=0.91399, beta =82.187, mu =-4.0011, lower.tail = TRUE,
  log.p = FALSE)


################################################################
#SPI-3 PROCEDURE#
################################################################
y=pre3
m=0
length(y)
n=length(y)
q=m/n  # m is total number of zero present in present time series
#n is total number of observations # q is probablity of zero.
q
Gx=CFit1
Gx
Hx=q+(1-q)*Gx
Hx
range(Hx)
#SPI Loop
sp=function(Hx){
  list=c()
  c0 = 2.515517
  c1 = 0.802583
  c2 = 0.010328
  d1 = 1.432788
  d2 = 0.189269
  d3 = 0.001308 
  for (i in 1:length(Hx)) {
    
    if((Hx[i]>0)&(Hx[i]<0.5)){
      
      t1=sqrt(log(((1/Hx[i])^2)))
      m1=((c0+c1*(t1)+c2*(t1^2))/(1+d1*(t1)+d2*(t1^2)+d3*(t1^3)))
      spiz=-(t1-m1)
      list[[i]]=spiz
      
    } else if((Hx[i]>0.5)&(Hx[i]<1.0)){
      
      t2=sqrt(log((1/(1-Hx[i])^2)))
      m2=((c0+c1*(t2)+c2*(t2^2))/(1+d1*(t2)+d2*(t2^2)+d3*(t2^3)))
      spiz=+(t2-m2)
      list[[i]]=spiz
    }
    
    
  } 
  
  print(list)
}
class(sp(Hx))
SPI3=sp(Hx)

plot.ts(SPI3)
spi3=SPI3
spi3
year=rep(c(1990:2018),each=12)
#For ensembling data:
temp=read.csv(file.choose())
TT=temp[,-1]
temp=as.matrix(TT)
Temp=as.numeric(temp)
Temp
DF1=data.frame(year,pre,Temp,pre3,spi3)
DF1
write.csv(DF1,file="Ensemble data Burkhan.csv")
#######################
#Burkhan Anlysis report:
#Improve the predictiblity of monthly drought Condition:


#1990-2018 series #348 observation:
data1=read.csv(file.choose())
str(data1)
#fix(data1)

##########################################################################
#Temperature Forecasting

T1=data1$Temp

#train data  1990-2011 #22 year series #264 obs
#test data  2012-2018 #7 year series #84 obs


Train.Temp=T1[1:264]  #training set
Test.Temp=T1[265:348] #testing set
#Combine Test and Train data
Com.Temp=c(Train.Temp,Test.Temp)

TTs=ts(Com.Temp,start=c(1990,1),frequency=12)
library(forecast)
ATM=auto.arima(TTs) #Auto temp model #ARIMA(2,0,1)(1,1,0)[12]
summary(ATM)

obs=Test.Temp

ARF<- rep(0,84)
for(i in 1:84){
  Tseries<- TTs[1:(264+i)]
  Tfit = arima(Tseries, order=c(1,0,0),seasonal=c(0,0,0),include.mean=TRUE)
  ARF[i] <- forecast(Tfit, h=1)$mean
}
summary(Tfit) #AIC 1329.38
Tem.Forecast=ARF   #forecast values
#Check Performance
MAPE <- mean(abs((ARF-obs)/obs))*100
RMSE<- sqrt(mean((ARF-obs)^2))
c(MAPE,RMSE)    

getPerformance = function(pred, val){
  n=length(val)
  res = pred - val
  MAE = sum(abs(res))/length(val)
  RSS = sum(res^2)
  MSE = RSS/length(val)
  RMSE = sqrt(MSE)
  MAPE=(1/n)*(sum(100 *(abs((val - pred) / val))))
  perf = data.frame(MAE, RSS, MSE, RMSE,MAPE)
  return(perf)
}              
getPerformance(ARF,obs)

#Resullt
#order=c(1,0,0),seasonal=c(0,0,0) #AIC 1329.38
#MAPE=0.4369448 #MSE=0.0132124  #RMSE=0.1149452 #MAE=0.1013614 

acf(Tfit$residuals,lag=2)
shapiro.test(Tfit$residuals)
#hist(Tfit$residuals)

T.Residuals=Tfit$residuals #Residuals from temperature model  #Whole series
Fitted.Temp=forecast(Tfit, h=1)$fitted #Whole series

############################################################################
#SPI3 Forecasting

str(data1)
#SS1=data1$spi3
#SS1=na.omit(SS1)
#SS1=c(0,0,SS1)

#train data  1990-2011 #22 year series #264 obs
#test data  2012-2018 #7 year series #84 obs

spi3=data1$spi3 
#spi3=na.omit(spi3) 
#spi3=c(0,0,spi3) 

Train.SP3=spi3[1:264] #training set
Test.SP3=spi3[265:348] #testing set
# n1= 240
# n2=72
Com.SP3=c(Train.SP3,Test.SP3)
TSs=ts(Com.SP3,start=c(1990,1),frequency=12)
#plot.ts(ser)
library(forecast)
M1=auto.arima(TSs)
summary(M1)
obs=Test.SP3
ARF<- rep(0,84)
for(i in 1:84){
  Sseries<- TSs[1:(264+i)]
  Sfit = arima(Sseries, order=c(0,2,1),seasonal=c(0,0,0),include.mean=TRUE)
  ARF[i] <- forecast(Sfit, h=1)$mean
}
summary(Sfit)
SPI3.Forecast=ARF
MAPE <- mean(abs((ARF-obs)/obs))*100
RMSE<- sqrt(mean((ARF-obs)^2))
c(MAPE,RMSE) #More best forecast Result

getPerformance = function(pred, val){
  n=length(val)
  res = pred - val
  MAE = sum(abs(res))/length(val)
  RSS = sum(res^2)
  MSE = RSS/length(val)
  RMSE = sqrt(MSE)
  MAPE=(1/n)*(sum(100 *(abs((val - pred) / val))))
  perf = data.frame(MAE, RSS, MSE, RMSE,MAPE)
  return(perf)
}              
getPerformance(ARF,obs)


#Resullt
#(order=c(0,2,1),seasonal=c(0,0,0)) #AIC 597.14
#MAPE=0.3650634 #MSE=1.340146e-05 #RMSE=0.003660801 #MAE=0.002933097 

acf(Tfit$residuals,lag=2)
shapiro.test(Tfit$residuals)
#hist(Tfit$residuals)

S.Residuals=Sfit$residuals #Residuals from temperature model  #Whole series
Fitted.Spi3=forecast(Sfit, h=1)$fitted #Whole series


############################################################################
# 3 Month moving averge Precipation Modeling

Pre3=data1$pre3

#train data  1990-2009 #20 year series #240 obs
#test data  2010-2015 #6 year series #72 obs
#Check normality of data

library("forecast")
#lambda=BoxCox.lambda(Pre3)
#lambda
#P3.Trans=BoxCox(Pre3,lambda)
#hist(P3.Trans)
#shapiro.test(trans)

#Using Log Transformation
Pre33=Pre3[3:348]
P3.Trans=log(Pre33)
#hist(P3.Trans)
P3.Trans=c(0,0,P3.Trans)

Train.Pre3=P3.Trans[1:264]  #training set 
Test.Pre3=P3.Trans[265:348] #testing set

#Combine Test and Train data
Com.Pre3=c(Train.Pre3,Test.Pre3)
#hist(Com.Pre3)

PTs=ts(Com.Pre3,start=c(1990,3),frequency=12)
#acf(PTs)
#pacf(PTs)
#dd=diff(PTs,2)
#acf(dd)
#pacf(dd)
#ddi=diff(PTs,difference=12)
#acf(ddi)
#pacf(ddi)

library(forecast)
APM=auto.arima(PTs) #ARIMA(5,0,2)(1,0,0)[12] with non-zero mean 
summary(APM)

obs=Test.Pre3
ARF<- rep(0,84)
for(i in 1:84){
  Pseries<- PTs[1:(264+i)]
  Pfit = arima(Pseries, order=c(0,0,0),seasonal=c(0,2,1),include.mean=TRUE)
  ARF[i] <- forecast(Pfit, h=1)$mean
}
summary(Pfit) #AIC 611.6
Pre3.Forecast=ARF   #forecast values
#Check Performance
MAPE <- mean(abs((ARF-obs)/obs))*100
RMSE<- sqrt(mean((ARF-obs)^2))
c(MAPE,RMSE)    

getPerformance = function(pred, val){
  n=length(val)
  res = pred - val
  MAE = sum(abs(res))/length(val)
  RSS = sum(res^2)
  MSE = RSS/length(val)
  RMSE = sqrt(MSE)
  MAPE=(1/n)*(sum(100 *(abs((val - pred) / val))))
  perf = data.frame(MAE, RSS, MSE, RMSE,MAPE)
  return(perf)
}              
getPerformance(ARF,obs)

#Resullt
#order=c(0,0,0),seasonal=c(0,2,1) #AIC 611.6
#MAPE=0.3650634 #RMSE=0.01637737 # MAE=0.01587033 #MSE=0.0002682182 

acf(Pfit$residuals,lag=2)
shapiro.test(Pfit$residuals)
#hist(Pfit$residuals) #Normaly Distributed

P.Residuals=Pfit$residuals #Residuals from Pre3 model  #Whole series
Fitted.Pre3=forecast(Pfit, h=1)$fitted  #Whole series

#############################################################################

#For climate werighing scheme (integrated climate predictor)

predicted.prec




#########################################################################

#Seclect fitted values of model only in the range of train data (1990-2009)
#For making of fit plolynoimail regression plot #240 observation


###########################################
#Temperature Fitted value
###########################################

Fitted.Temp=forecast(Tfit, h=1)$fitted #Whole series

fit1.temp=ts(Fitted.Temp,start=c(1990,1),end=c(2011,12),frequency=12)
#Get all data from class of "ts"
#current Month is "April"
fitt.april=subset(fit1.temp, cycle(fit1.temp) == 4)
fitt.april
#Target Months
fitt.may=subset(fit1.temp, cycle(fit1.temp) == 5)
fitt.may
fitt.jun=subset(fit1.temp, cycle(fit1.temp) == 6)
fitt.jun
fitt.jul=subset(fit1.temp, cycle(fit1.temp) == 7)
fitt.jul
par(mfrow=c(2,2))
library("stats")
#scatter plot with locally weighted polynoimail regression
#scatter.smooth(fitt.april,fitt.may,xlab = "Mar",ylab="May",main="Temperature/lag-1") 
cor(fitt.april,fitt.may) #0.9707566
#scatter.smooth(fitt.april,fitt.jun,xlab = "Mar",ylab="Jun",main="Temperature/lag-2") 
cor(fitt.april,fitt.jun)  #0.9760443
#scatter.smooth(fitt.april,fitt.jul,xlab = "Mar",ylab="Jul",main="Temperature/lag-3") 
cor(fitt.april,fitt.jul)   #0.961564

#########################################
#pre3 fitted values
#########################################

Fitted.Pre3=forecast(Pfit, h=1)$fitted  #Whole series

pre.ts=ts(Fitted.Pre3,start=c(1990,1),end=c(2011,12),frequency=12)

fitp.april=subset(pre.ts, cycle(pre.ts) == 4)
fitp.april
#class(fitp.april)
fitp.may=subset(pre.ts, cycle(pre.ts) == 5)
fitp.may
fitp.jun=subset(pre.ts, cycle(pre.ts) == 6)
fitp.jun
fitp.jul=subset(pre.ts, cycle(pre.ts) == 7)
fitp.jul
par(mfrow=c(2,2))
library("stats")
#scatter plot with locally weighted polynoimail regression
scatter.smooth(fitp.april,fitp.may,xlab = "Mar",ylab="May",main="Pre-3/lag-1") 
cor(fitp.april,fitp.may) #0.7345049
scatter.smooth(fitp.april,fitp.jun,xlab = "Mar",ylab="Jun",main="Pre-3/lag-2") 
cor(fitp.april,fitp.jun)  #0.6436559
scatter.smooth(fitp.april,fitp.jul,xlab = "Mar",ylab="Jul",main="Pre-3/lag-3") 
cor(fitp.april,fitp.jul)  #0.5969224

######################################################
#Spi3 fitted Values
######################################################
Fitted.Spi3=forecast(Sfit, h=1)$fitted #Whole series

fit1.spi=ts(Fitted.Spi3,start=c(1990,1),end=c(2011,12),frequency=12)

#Get all data from class of "ts"

#current Month is "March"
fits.april=subset(fit1.spi, cycle(fit1.spi) == 4)
fits.april
#Target Months
fits.may=subset(fit1.spi, cycle(fit1.spi) == 5)
fits.may
fits.jun=subset(fit1.spi, cycle(fit1.spi) == 6)
fits.jun
fits.jul=subset(fit1.spi, cycle(fit1.spi) == 7)
fits.jul
par(mfrow=c(2,2))
library("stats")
#scatter plot with locally weighted polynoimail regression
scatter.smooth(fits.april,fits.may,xlab = "Mar",ylab="May",main="Spi/lag-1") 
cor(fits.april,fits.may) #0.5969224
scatter.smooth(fits.april,fits.jun,xlab = "Mar",ylab="May",main="Spi/lag-2") 
cor(fits.april,fits.jun)   #0.6618973
scatter.smooth(fits.april,fits.jul,xlab = "Mar",ylab="Jun",main="Spi/lag-3") 
cor(fits.april,fits.jul)    #0.5700374

################################################################################
         #Fitted value of Model ggplot
################################################################################


library("ggplot2")
library("dplyr")
#FOR SPI Lag Plots
dats = data.frame(fits.april,fits.may,fits.jun,fits.jul)
cor(fits.april,fits.may)   #0.7532025
range(fits.may)
s1=ggplot(dats, aes(fits.april, fits.may)) + geom_point(size=4, shape=20)+
  geom_smooth(method=loess, se=FALSE, fullrange=TRUE)+
  scale_y_continuous(limits = c(-3,3))+
  geom_text(x=-2.3, y=2.1, label="r=0.75", color="black",fontface="italic",hjust=0.5)+
  labs(x ="April", y = "May", size=12, face="bold.italic")+
  theme(
    plot.title = element_text(color="black", size=12, face="bold.italic"),
    axis.title.x = element_text(color="black", size=10, face="plain"),
    axis.title.y = element_text(color="black", size=10, face="plain"))+
  facet_grid(. ~ "SPI-3 /Lag-1")+
  theme(strip.text.x = element_text(face="plain", size=12, colour="black")) +
  theme(strip.text.y = element_text(face="plain", size=12, color="white")) +
  theme(strip.background = element_rect(fill="grey"))    
s1

cor(fits.april,fits.jun)  #0.6618973
s2=ggplot(dats, aes(fits.april, fits.jun)) + geom_point(size=4, shape=20)+
  geom_smooth(method=loess, se=FALSE, fullrange=TRUE)+
scale_y_continuous(limits = c(-3,3))+
  geom_text(x=-2.3, y=2.1, label="r=0.66", color="black",fontface="italic",hjust=0.5)+
  labs(x ="April", y = "Jun", size=12, face="bold.italic")+
  theme(
    plot.title = element_text(color="black", size=12, face="bold.italic"),
    axis.title.x = element_text(color="black", size=10, face="plain"),
    axis.title.y = element_text(color="black", size=10, face="plain"))+
    facet_grid(. ~ "SPI-3 /Lag-2")+
  theme(strip.text.x = element_text(face="plain", size=12, colour="black")) +
  theme(strip.text.y = element_text(face="plain", size=12, color="white")) +
  theme(strip.background = element_rect(fill="grey"))
s2

cor(fits.april,fits.jul) #0.5700374
s3=ggplot(dats, aes(fits.april, fits.jul)) + geom_point(size=4, shape=20)+
  geom_smooth(method=loess, se=FALSE, fullrange=TRUE)+
scale_y_continuous(limits = c(-3,3))+
  geom_text(x=-2.3, y=2.1, label="r=0.57", color="black",fontface="italic",hjust=0.5)+
  labs(x ="April", y = "July", size=12, face="bold.italic")+
  theme(
    plot.title = element_text(color="black", size=12, face="bold.italic"),
    axis.title.x = element_text(color="black", size=10, face="plain"),
    axis.title.y = element_text(color="black", size=10, face="plain"))+
    facet_grid(. ~ "SPI-3 /Lag-3")+
  theme(strip.text.x = element_text(face="plain", size=12, colour="black")) +
  theme(strip.text.y = element_text(face="plain", size=12, color="white")) +
  theme(strip.background = element_rect(fill="grey"))
s3

#Combine spi ggplots into one page
library(ggpubr)
figure1 <- ggarrange(s1, s2, s3, ncol = 2, nrow = 2)
figure1                    
                   
#For Temperature
datt = data.frame(fitt.april,fitt.may,fitt.jun,fitt.jul)
cor(fitt.april,fitt.may)  #0.9707566
range(fitt.may)
t1=ggplot(datt, aes(fitt.april, fitt.may)) + geom_point(size=4, shape=20)+
  geom_smooth(method=loess, se=FALSE, fullrange=TRUE)+
scale_y_continuous(limits = c(10,36))+
  geom_text(x=13.7, y=32.9, label="r=0.97", color="black",fontface="italic",hjust=0.5)+
  labs(x ="April", y = "May", size=12, face="bold.italic")+
  theme(
    plot.title = element_text(color="black", size=12, face="bold.italic"),
    axis.title.x = element_text(color="black", size=10, face="plain"),
    axis.title.y = element_text(color="black", size=10, face="plain"))+
    facet_grid(. ~ "Temperature /Lag-1")+
  theme(strip.text.x = element_text(face="plain", size=12, colour="black")) +
  theme(strip.text.y = element_text(face="plain", size=12, color="white")) +
  theme(strip.background = element_rect(fill="grey"))
t1


cor(fitt.april,fitt.jun)   #0.9760443
range(fitt.jun)
t2=ggplot(datt, aes(fitt.april, fitt.jun)) + geom_point(size=4, shape=20)+
  geom_smooth(method=loess, se=FALSE, fullrange=TRUE)+
scale_y_continuous(limits = c(10,36))+
  geom_text(x=13.7, y=32.9, label="r=0.97", color="black",fontface="italic",hjust=0.5)+
  labs(x ="April", y = "Jun", size=12, face="bold.italic")+
  theme(
    plot.title = element_text(color="black", size=12, face="bold.italic"),
    axis.title.x = element_text(color="black", size=10, face="plain"),
    axis.title.y = element_text(color="black", size=10, face="plain"))+
    facet_grid(. ~ "Temperature /Lag-2")+
  theme(strip.text.x = element_text(face="plain", size=12, colour="black")) +
  theme(strip.text.y = element_text(face="plain", size=12, color="white")) +
  theme(strip.background = element_rect(fill="grey"))
t2

cor(fitt.april,fitt.jul)   #0.961564
range(fitt.jul)
t3=ggplot(datt, aes(fitt.april, fitt.jul)) + geom_point(size=4, shape=20)+
  geom_smooth(method=loess, se=FALSE, fullrange=TRUE)+
scale_y_continuous(limits = c(10,36))+
  geom_text(x=13.7, y=32.9, label="r=0.96", color="black",fontface="italic",hjust=0.5)+
  labs(x ="April", y = "July", size=12, face="bold.italic")+
  theme(
    plot.title = element_text(color="black", size=12, face="bold.italic"),
    axis.title.x = element_text(color="black", size=10, face="plain"),
    axis.title.y = element_text(color="black", size=10, face="plain"))+
    facet_grid(. ~ "Temperature /Lag-3")+
  theme(strip.text.x = element_text(face="plain", size=12, colour="black")) +
  theme(strip.text.y = element_text(face="plain", size=12, color="white")) +
  theme(strip.background = element_rect(fill="grey"))
t3

#Combine temperature ggplots into one page
library(ggpubr)
figure2 <- ggarrange(t1, t2, t3, ncol = 2, nrow = 2)
figure2
#For Preciaption fitted plots
datp = data.frame(fitp.april,fitp.may,fitp.jun,fitp.jul)
cor(fitp.april,fitp.may)   #0.7345049
range(fitp.may)
p1=ggplot(datp, aes(fitp.april, fitp.may)) + geom_point(size=4, shape=20)+
  geom_smooth(method=loess, se=FALSE, fullrange=TRUE)+
scale_y_continuous(limits = c(1,8))+
  geom_text(x=2.3, y=7.1, label="r=0.73", color="black",fontface="italic",hjust=0.5)+
  labs(x ="April", y = "May", size=12, face="bold.italic")+
  theme(
    plot.title = element_text(color="black", size=12, face="bold.italic"),
    axis.title.x = element_text(color="black", size=10, face="plain"),
    axis.title.y = element_text(color="black", size=10, face="plain"))+
    facet_grid(. ~ "Precipitation-3 /Lag-1")+
  theme(strip.text.x = element_text(face="plain", size=12, colour="black")) +
  theme(strip.text.y = element_text(face="plain", size=12, color="white")) +
  theme(strip.background = element_rect(fill="grey"))
p1

cor(fitp.april,fitp.jun)    # 0.6436559
range(fitp.jun)
p2=ggplot(datp, aes(fitp.april, fitp.jun)) + geom_point(size=4, shape=20)+
  geom_smooth(method=loess, se=FALSE, fullrange=TRUE)+
scale_y_continuous(limits = c(1,8))+
  geom_text(x=2.3, y=7.1, label="r=0.64", color="black",fontface="italic",hjust=0.5)+
  labs(x ="April", y = "Jun", size=12, face="bold.italic")+
  theme(
    plot.title = element_text(color="black", size=12, face="bold.italic"),
    axis.title.x = element_text(color="black", size=10, face="plain"),
    axis.title.y = element_text(color="black", size=10, face="plain"))+
    facet_grid(. ~ "Precipitation-3 /Lag-2")+
  theme(strip.text.x = element_text(face="plain", size=12, colour="black")) +
  theme(strip.text.y = element_text(face="plain", size=12, color="white")) +
  theme(strip.background = element_rect(fill="grey"))
p2

cor(fitp.april,fitp.jul)  # 0.5969224
range(fitp.jul)
p3=ggplot(datp, aes(fitp.april, fitp.jul)) + geom_point(size=4, shape=20)+
  geom_smooth(method=loess, se=FALSE, fullrange=TRUE)+
 scale_y_continuous(limits = c(1,8))+
  geom_text(x=2.3, y=7.1, label="r= 0.59", color="black",fontface="italic",hjust=0.5)+
  labs(x ="April", y = "July", size=12, face="bold.italic")+ 
  theme(
    plot.title = element_text(color="black", size=12, face="bold.italic"),
    axis.title.x = element_text(color="black", size=10, face="plain"),
    axis.title.y = element_text(color="black", size=10, face="plain"))+
    facet_grid(. ~ "Precipitation-3 /Lag-3")+
  theme(strip.text.x = element_text(face="plain", size=12, colour="black")) +
  theme(strip.text.y = element_text(face="plain", size=12, color="white")) +
  theme(strip.background = element_rect(fill="grey"))
p3

#Combine precipation ggplots into one page
library(ggpubr)
figure3 <- ggarrange(p1, p2, p3, ncol = 2, nrow = 2)
figure3
#Combine all ggplots into one page
library(ggpubr)
figure <- ggarrange(s1,t1,p1,s2,t2,p2,s3, t3, p3, ncol = 3, nrow = 3)
figure

###########################################################################################
#Oh this stage all work done on the Test data beacause comparison of obs and predicted:

#######################################
# Forecasted test data 2010-2015 (72 obs)
#######################################
#Convert SPI into time series data
SPI3.Forecast
#Convert into monthly data and select current and target month
#Convert into 12 month data frame
set.seed(12)
#check the dimemsion of data no of rows and colume 

#we convert into data frame
sp3.df <- as.data.frame(matrix(SPI3.Forecast,nrow=7,ncol=12))
sp3.df
rownames(sp3.df) <- seq(from=2012, to=2018)
colnames(sp3.df) <- c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")
sp3.df.m <- ts(as.vector(as.matrix(sp3.df)), start=c(2012,1), end=c(2018,12), frequency=12)
sp3.df.m 
#seclet data
jan.sp3=subset(sp3.df.m,cycle(sp3.df.m)==1)
feb.sp3=subset(sp3.df.m,cycle(sp3.df.m)==2)
mar.sp3=subset(sp3.df.m,cycle(sp3.df.m)==3)
apr.sp3=subset(sp3.df.m,cycle(sp3.df.m)==4)
may.sp3=subset(sp3.df.m,cycle(sp3.df.m)==5)
jun.sp3=subset(sp3.df.m,cycle(sp3.df.m)==6)
jul.sp3=subset(sp3.df.m,cycle(sp3.df.m)==7)
aug.sp3=subset(sp3.df.m,cycle(sp3.df.m)==8)
sep.sp3=subset(sp3.df.m,cycle(sp3.df.m)==9)
oct.sp3=subset(sp3.df.m,cycle(sp3.df.m)==10)
nov.sp3=subset(sp3.df.m,cycle(sp3.df.m)==11)
dec.sp3=subset(sp3.df.m,cycle(sp3.df.m)==12)
#we take current month "April"
apr.sp3=subset(sp3.df.m,cycle(sp3.df.m)==4)

########################################################
          #For taget month "May"
########################################################
#Fit local regression model
library("locfit")
may.lgr=locfit(may.sp3~apr.sp3,alpha=0.7)
summary(may.lgr)
res.lgr.may=residuals(may.lgr) #Residuals from locfit model
#fit cubic hermitte 
cubic.hermitte.may=predict(may.lgr, newdata=NULL, where="data", tr=NULL, what="infl",
            band="none", get.data=TRUE,se.fit=TRUE)
cubic.hermitte.may
#plot.ts(cubic.hermitte.may$fit)
range(cubic.hermitte.may$fit)
#y00 is intial mean forecast of May target month:
y00.may=cubic.hermitte.may$fit
y00.may

###################################################
#save locfit model residules and resmapling "May"
###################################################
res.lgr.may  #residuals from locfit model target month may
year=rep(c(2012:2018), each = 1)
res.df.may=data.frame(year,res.lgr.may)
res.df.may
#Resampling through Bootstraps
library("rsample")
library("tidyr")
samples = 1000 #times
boots = bootstraps(res.df.may, times=samples)
head(boots)
#split into 1000 times
bootstrap_sample = analysis(boots$splits[[1]]) #first split
mean(bootstrap_sample[,2])  #mean of first split
bootstrap_sample
head(bootstrap_sample)
mean(analysis(boots$splits[[1]])$res.lgr.may) #first split mean
bootstrap_sample = analysis(boots$splits[[8]]) #select 8th sample 
#hist(bootstrap_sample$res.lgr.apr)
#hist(res.df.apr$res.lgr.apr)
bootstrap_sample
ee1.may=bootstrap_sample$res.lgr.may  #Resample residuals from locfit model
ee1.may #Resample residuals from locfit model apr


########################################################
          #For taget month "Jun"
########################################################
#Fit local regression model
library("locfit")
jun.lgr=locfit(jun.sp3~apr.sp3,alpha=0.7)
summary(jun.lgr)
res.lgr.jun=residuals(jun.lgr) #Residuals from locfit model
#fit cubic hermitte 
cubic.hermitte.jun=predict(jun.lgr, newdata=NULL, where="data", tr=NULL, what="infl",
            band="none", get.data=TRUE,se.fit=TRUE)
cubic.hermitte.jun
#plot.ts(cubic.hermitte.jun$fit)
range(cubic.hermitte.jun$fit)
#y00 is intial mean forecast of target month:
y00.jun=cubic.hermitte.jun$fit
y00.jun
###################################################
#save locfit model residules and resmapling "Jun"
###################################################
res.lgr.jun  #residuals from locfit model target month jun
year=rep(c(2012:2018), each = 1)
res.df.jun=data.frame(year,res.lgr.jun)
res.df.jun
#Resampling through Bootstraps
library("rsample")
library("tidyr")
samples = 1000 #times
boots = bootstraps(res.df.jun, times=samples)
head(boots)
#split into 1000 times
bootstrap_sample = analysis(boots$splits[[1]]) #first split
mean(bootstrap_sample[,2])  #mean of first split
bootstrap_sample
head(bootstrap_sample)
mean(analysis(boots$splits[[1]])$res.lgr.jun) #first split mean
bootstrap_sample = analysis(boots$splits[[8]]) #select 8th sample 
#hist(bootstrap_sample$res.lgr.jun)
#hist(res.df.may$res.lgr.jun)
bootstrap_sample
ee1.jun=bootstrap_sample$res.lgr.jun  #Resample residuals from locfit model
ee1.jun #Resample residuals from locfit model jun

########################################################
          #For taget month "Jul"
########################################################
#Fit local regression model
library("locfit")
jul.lgr=locfit(jul.sp3~apr.sp3,alpha=0.7)
summary(jul.lgr)
res.lgr.jul=residuals(jul.lgr) #Residuals from locfit model
#fit cubic hermitte 
cubic.hermitte.jul=predict(jul.lgr, newdata=NULL, where="data", tr=NULL, what="infl",
            band="none", get.data=TRUE,se.fit=TRUE)
cubic.hermitte.jul
#plot.ts(cubic.hermitte.jul$fit)
range(cubic.hermitte.jul$fit)
#y00 is intial mean forecast of target month jul:
y00.jul=cubic.hermitte.jul$fit
y00.jul
###################################################
#save locfit model residules and resmapling "July"
###################################################
res.lgr.jul  #residuals from locfit model target month jun
year=rep(c(2012:2018), each = 1)
res.df.jul=data.frame(year,res.lgr.jul)
res.df.jul
#Resampling through Bootstraps
library("rsample")
library("tidyr")
samples = 1000 #times
boots = bootstraps(res.df.jul, times=samples)
head(boots)
#split into 1000 times
bootstrap_sample = analysis(boots$splits[[1]]) #first split
mean(bootstrap_sample[,2])  #mean of first split
bootstrap_sample
head(bootstrap_sample)
mean(analysis(boots$splits[[1]])$res.lgr.jul) #first split mean
bootstrap_sample = analysis(boots$splits[[8]]) #select 8th sample 
#hist(bootstrap_sample$res.lgr.jul)
#hist(res.df.jun$res.lgr.jul)
bootstrap_sample
ee1.jul=bootstrap_sample$res.lgr.jul  #Resample residuals from locfit model
ee1.jul #Resample residuals from locfit model jun

##########################################################################################

T1
Tr1=T1[265:348]   #Select only Original Pre3 data in the range of test data values 
Tr1
#covert into 12 month data frame and select target month:
set.seed(12)
#check the dimemsion of data no of rows and colume 

#we convert into data frame
Tr1.df <- as.data.frame(matrix(Tr1,nrow=7,ncol=12))
Tr1.df
rownames(Tr1.df) <- seq(from=2012, to=2018)
colnames(Tr1.df) <- c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")
tem.df.m <- ts(as.vector(as.matrix(Tr1.df)), start=c(2012,1), end=c(2018,12), frequency=12)
tem.df.m
#select month
jan.tem=subset(tem.df.m,cycle(tem.df.m)==1)
feb.tem=subset(tem.df.m,cycle(tem.df.m)==2)
mar.tem=subset(tem.df.m,cycle(tem.df.m)==3)
apr.tem=subset(tem.df.m,cycle(tem.df.m)==4)
may.tem=subset(tem.df.m,cycle(tem.df.m)==5)
jun.tem=subset(tem.df.m,cycle(tem.df.m)==6)
jul.tem=subset(tem.df.m,cycle(tem.df.m)==7)
aug.tem=subset(tem.df.m,cycle(tem.df.m)==8)
sep.tem=subset(tem.df.m,cycle(tem.df.m)==9)
oct.tem=subset(tem.df.m,cycle(tem.df.m)==10)
nov.tem=subset(tem.df.m,cycle(tem.df.m)==11)
dec.tem=subset(tem.df.m,cycle(tem.df.m)==12)

Pre3
pr1=Pre3[265:348]   #Select only Original Pre3 data in the range of test data values 
pr1
#covert into 12 month data frame and select target month:
set.seed(12)
#check the dimemsion of data no of rows and colume 

#we convert into data frame
pr1.df <- as.data.frame(matrix(pr1,nrow=7,ncol=12))
pr1.df
rownames(pr1.df) <- seq(from=2012, to=2018)
colnames(pr1.df) <- c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")
pr3.df.m <- ts(as.vector(as.matrix(pr1.df)), start=c(2012,1), end=c(2018,12), frequency=12)
pr3.df.m 
#select month
jan.pre3=subset(pr3.df.m,cycle(pr3.df.m)==1)
feb.pre3=subset(pr3.df.m,cycle(pr3.df.m)==2)
mar.pre3=subset(pr3.df.m,cycle(pr3.df.m)==3)
apr.pre3=subset(pr3.df.m,cycle(pr3.df.m)==4)
may.pre3=subset(pr3.df.m,cycle(pr3.df.m)==5)
jun.pre3=subset(pr3.df.m,cycle(pr3.df.m)==6)
jul.pre3=subset(pr3.df.m,cycle(pr3.df.m)==7)
aug.pre3=subset(pr3.df.m,cycle(pr3.df.m)==8)
sep.pre3=subset(pr3.df.m,cycle(pr3.df.m)==9)
oct.pre3=subset(pr3.df.m,cycle(pr3.df.m)==10)
nov.pre3=subset(pr3.df.m,cycle(pr3.df.m)==11)
dec.pre3=subset(pr3.df.m,cycle(pr3.df.m)==12)



##########################################################################
                       #Work For May
##########################################################################
may.tem
length(may.tem)
ss=sort(may.tem)
ss
rank(ss)
r1=range(ss[1:2])
r1
r2=range(ss[3:4])
r2
r3=range(ss[5:7])
r3
w=function(x){
  
  df=NULL    
  for (i in 1:length(x)){
    if((x[i]>17.34)&(x[i]<17.71)){
      y=x[i]
      fac=c("Below")
      df = rbind(df,data.frame(y,fac))
    }else if((x[i]>17.70)&(x[i]<20.11)){
      y=x[i]
      fac=c("Normal")
      df = rbind(df,data.frame(y,fac))
      
    }else if((x[i]>20.10)&(x[i]<78.5)){
      
      y=x[i]
      fac=c("Above")
      df = rbind(df,data.frame(y,fac))
      
    }
    
  }
  return(df)
}

output.tem.may=w(may.tem)  #Categorize into level
output.tem.may

#########################################################
#Ordering data with respect to factor TEMPERATURE "May"
#########################################################
T.Residuals=Tfit$residuals #Residuals from temperature model  #Whole series
T.Residuals    #temperature model Residuals
length(T.Residuals)
tt.residual=ts(T.Residuals[265:348],start=c(2012,1),end=c(2018,12),frequency=12)
#tt.residual=ts(T.Residuals,start=c(1990,1),end=c(2015,12),frequency=12)
res.tem.may=subset(tt.residual,cycle(tt.residual)==5)
res.tem.may
year=rep(c(2012:2018), each = 1)
app=data.frame(year,output.tem.may,res.tem.may)
df.tem.may<- app[order(app$fac),]
df.tem.may #Order wise data

##########################################################
#Conditional Residuls resampling of temperature "May
##########################################################
df.tem.may
temp.resample=df.tem.may[,-2]
#For "Below Temperature" 
t.blow.may=temp.resample[6:7,]
mean(t.blow.may$res.tem.may)
library(rsample)
library(tidyverse)
samples = 300 #times
boots.tb.may = bootstraps(t.blow.may, times=samples)
head(boots.tb.may)
#split into 300 times
bootstrap_sample.tb.may = analysis(boots.tb.may$splits[[1]]) #first split
mean(bootstrap_sample.tb.may[,3])  #mean of first split
head(bootstrap_sample.tb.may)
mean(analysis(boots.tb.may $splits[[1]])$res.tem.may) #first split mean
#all split means in the colume of mn
boots_mean.tb.may = boots.tb.may%>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.tem.may)))
boots_mean.tb.may
#Approximated distribution of mean
#boots_mean.tb.may%>%
  ggplot(aes(mn)) +
  geom_density() +
  labs(title = "Approximated sampling distribution of mean res.tem.may",
       x = "mean res.tem.may",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.tb.may = analysis(boots.tb.may$splits[[12]]) #select 12 split or sample
#hist(bootstrap_sample.tb.may$res.tem.may)
#hist(t.blow.may$res.tem.may)
#For "Normal" Temperature "April"
temp.resample=df.tem.may[,-2]
temp.resample
t.normal.may=temp.resample[1:2,]
mean(t.normal.may$res.tem.may)
library(rsample)
library(tidyverse)
samples = 300 #times
boots.tn.may= bootstraps(t.normal.may, times=samples)
head(boots.tn.may)
#split into 300 times
bootstrap_sample.tn.may = analysis(boots.tn.may$splits[[1]]) #first split
mean(bootstrap_sample.tn.may[,3])  #mean of first split
head(bootstrap_sample.tn.may)
mean(analysis(boots.tn.may$splits[[1]])$res.tem.may) #first split mean
#all split means in the colume of mn
boots_mean.tn.may=boots.tn.may %>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.tem.may)))
boots_mean.tn.may
#Approximated distribution of mean
#boots_mean.tn.may%>%
  ggplot(aes(mn)) +
  geom_density() +
  labs(title = "Approximated sampling distribution of mean res.tem.may",
       x = "mean res.tem.may",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.tn.may= analysis(boots.tn.may$splits[[197]]) #select 197 split or sample
#hist(bootstrap_sample.tn.may$res.tem.may)
#hist(t.normal.may$res.tem.may)
#For "Above" Temperature
temp.resample
t.above.may=temp.resample[3:5,]
t.above.may
mean(t.above.may$res.tem.may)
library(rsample)
library(tidyverse)
samples = 300 #times
boots.ta.may= bootstraps(t.above.may, times=samples)
head(boots.ta.may)
#split into 300 times
bootstrap_sample.ta.may = analysis(boots.ta.may$splits[[1]]) #first split
mean(bootstrap_sample.ta.may[,3])  #mean of first split
head(bootstrap_sample.ta.may)
mean(analysis(boots.ta.may$splits[[1]])$res.tem.may) #first split mean
#all split means in the colume of mn
boots_mean.ta.may=boots.ta.may %>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.tem.may)))
boots_mean.ta.may
#Approximated distribution of mean
#boots_mean.ta.may%>%
  ggplot(aes(mn)) +
  geom_density() +
 labs(title = "Approximated sampling distribution of mean res.tem.may",
       x = "mean res.tem.may",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.ta.may= analysis(boots.ta.may$splits[[156]]) #select 156 split or sample
#hist(bootstrap_sample.ta.may$res.tem.may)
#hist(t.above.may$res.tem.may)
#Combine all categories of residuals into one data frame:
t.residuals.may=rbind(bootstrap_sample.tb.may,bootstrap_sample.tn.may,bootstrap_sample.ta.may)
mean(t.residuals.may$res.tem.may)
#mean(temp$res.tem.may)
t.residuals.may #resample residuls

##################################################
#Precipation-3 Ranking target month May
##################################################
may.pre3
#Ranking
ss=sort(may.pre3)
ss
rank(ss)
r1=range(ss[1:2])
r1
r2=range(ss[3:4])
r2
r3=range(ss[5:7])
r3

w=function(x){
  
  df=NULL    
  for (i in 1:length(x)){
    if((x[i]>21.9)&(x[i]<29.80)){
      y=x[i]
      fac=c("Below")
      df = rbind(df,data.frame(y,fac))
    }else if((x[i]>29.79)&(x[i]<75.92)){
      y=x[i]
      fac=c("Normal")
      df = rbind(df,data.frame(y,fac))
      
    }else if((x[i]>75.91)){
      
      y=x[i]
      fac=c("Above")
      df = rbind(df,data.frame(y,fac))
      
    }
    
  }
  return(df)
}

output.pre3.may=w(may.pre3) #Categorize pre3 into levels
output.pre3.may
############################################################
#Ordering data with respect to factor Precipations-3 "May"
############################################################
P.Residuals=Pfit$residuals #Residuals from Pre3 model  #Whole series

pp.residual=ts(P.Residuals,start=c(2012,1),end=c(2018,12),frequency=12)

res.pre.may=subset(pp.residual,cycle(pp.residual)==5)
res.pre.may
year=rep(c(2012:2018), each = 1)
app=data.frame(year,output.pre3.may,res.pre.may)
df.pre3.may<- app[order(app$fac),]
df.pre3.may #Order wise data
##########################################################
#Conditional Residuls resampling of Precipation-3 "May
##########################################################
df.pre3.may #Order wise data

pre.3=df.pre3.may[,-2]
pre.3
#For "Below" Precipation3 
p.blow.may=pre.3[4:5,]
mean(p.blow.may$res.pre.may)
library(rsample)
library(tidyverse)
samples = 500 #times
boots.pb.may = bootstraps(p.blow.may, times=samples)
head(boots.pb.may)
#split into 300 times
bootstrap_sample.pb.may = analysis(boots.pb.may$splits[[1]]) #first split
mean(bootstrap_sample.pb.may[,3])  #mean of first split
head(bootstrap_sample.pb.may)
mean(analysis(boots.pb.may $splits[[1]])$res.pre.may) #first split mean
#all split means in the colume of mn
boots_mean.pb.may = boots.pb.may%>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.pre.may)))
boots_mean.pb.may
#Approximated distribution of mean
#boots_mean.pb.may%>%
  ggplot(aes(mn)) +
  geom_density() +
  labs(title = "Approximated sampling distribution of mean res.pre.may",
       x = "mean res.pre.may",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.pb.may = analysis(boots.pb.may$splits[[27]]) #select 27 split or sample
#hist(bootstrap_sample.pb.may$res.pre.may)
#hist(p.blow.may$res.pre.may)

#For "Normal" Preciaption3 
pre.3
p.normal.may=pre.3[6:7,]
mean(p.normal.may$res.pre.may)
library(rsample)
library(tidyverse)
samples = 500 #times
boots.pn.may = bootstraps(p.normal.may, times=samples)
head(boots.pn.may)
#split into 300 times
bootstrap_sample.pn.may = analysis(boots.pn.may$splits[[1]]) #first split
mean(bootstrap_sample.pn.may[,3])  #mean of first split
head(bootstrap_sample.pn.may)
mean(analysis(boots.pn.may $splits[[1]])$res.pre.may) #first split mean
#all split means in the colume of mn
boots_mean.pn.may = boots.pn.may%>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.may.apr)))
boots_mean.pn.may
#Approximated distribution of mean
#boots_mean.pn.may%>%
  ggplot(aes(mn)) +
  geom_density() +
  labs(title = "Approximated sampling distribution of mean res.pre.may",
       x = "mean res.pre.may",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.pn.may = analysis(boots.pn.may$splits[[29]]) #select 29 split or sample
#hist(bootstrap_sample.pn.may$res.pre.may)
#hist(p.normal.may$res.pre.may)

#For "Above" Preciaption3 
pre.3
p.above.may=pre.3[1:3,]
mean(p.normal.may$res.pre.may)
library(rsample)
library(tidyverse)
samples = 500 #times
boots.pa.may = bootstraps(p.above.may, times=samples)
head(boots.pa.may)
#split into 500 times
bootstrap_sample.pa.may = analysis(boots.pa.may$splits[[1]]) #first split
mean(bootstrap_sample.pa.may[,3])  #mean of first split
head(bootstrap_sample.pa.may)
mean(analysis(boots.pa.may $splits[[1]])$res.pre.may) #first split mean
#all split means in the colume of mn
boots_mean.pa.may = boots.pa.may%>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.pre.may)))
boots_mean.pa.may
#Approximated distribution of mean
#boots_mean.pa.may%>%
  ggplot(aes(mn)) +
  geom_density() +
  labs(title = "Approximated sampling distribution of mean res.pre.may",
       x = "mean res.pre.may",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.pa.may = analysis(boots.pa.may$splits[[77]]) #select 77 split or sample
#hist(bootstrap_sample.pa.may$res.pre.may)
#hist(p.above.may$res.pre.may)
#Combine all categories of residuals into one data frame:
p.residuals.may=rbind(bootstrap_sample.pb.may,bootstrap_sample.pa.may,bootstrap_sample.pn.may)
mean(p.residuals.may$res.pre.may)
mean(pre.3$res.pre.may)
p.residuals.may  #resample residuls of pre3 of May

########################################################################
             #End Work For May
########################################################################


#########################################################################
                        #Work for Jun
#########################################################################
#Rank the observed of target month "Jun" temperature and pre3 into three categories 
jun.tem

length(jun.tem)
ss=sort(jun.tem)
ss
rank(ss)
r1=range(ss[1:2])
r1
r2=range(ss[3:4])
r2
r3=range(ss[5:7])
r3
w=function(x){
  
  df=NULL    
  for (i in 1:length(x)){
    if((x[i]>16.1)&(x[i]<17.1)){
      y=x[i]
      fac=c("Below")
      df = rbind(df,data.frame(y,fac))
    }else if((x[i]>17.0)&(x[i]<23.56)){
      y=x[i]
      fac=c("Normal")
      df = rbind(df,data.frame(y,fac))
      
    }else if((x[i]>23.55)){
      
      y=x[i]
      fac=c("Above")
      df = rbind(df,data.frame(y,fac))
      
    }
    
  }
  return(df)
}

output.tem.jun=w(jun.tem)   #Categorize into level
output.tem.jun

#########################################################
#Ordering data with respect to factor TEMPERATURE "Jun"
#########################################################
T.Residuals=Tfit$residuals #Residuals from temperature model  #Whole series
length(T.Residuals)
tt.residual=ts(T.Residuals[265:348],start=c(2012,1),end=c(2018,12),frequency=12)
#tt.residual=ts(T.Residuals,start=c(1990,1),end=c(2015,12),frequency=12)
res.tem.jun=subset(tt.residual,cycle(tt.residual)==6)
res.tem.jun
year=rep(c(2012:2018), each = 1)
app=data.frame(year,output.tem.jun,res.tem.jun)
df.tem.jun<- app[order(app$fac),]
df.tem.jun#Order wise data

##########################################################
#Conditional Residuls resampling of temperature "Jun"
##########################################################
df.tem.jun
temp.jun=df.tem.jun[,-2]
#For "Below Temperature" 
t.blow.jun=temp.jun[6:7,]
mean(t.blow.jun$res.tem.jun)
library(rsample)
library(tidyverse)
samples = 300 #times
boots.tb.jun= bootstraps(t.blow.jun, times=samples)
head(boots.tb.jun)
#split into 300 times
bootstrap_sample.tb.jun= analysis(boots.tb.jun$splits[[1]]) #first split
mean(bootstrap_sample.tb.jun[,3])  #mean of first split
head(bootstrap_sample.tb.jun)
mean(analysis(boots.tb.jun$splits[[1]])$res.tem.jun) #first split mean
#all split means in the colume of mn
boots_mean.tb.jun= boots.tb.jun%>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.tem.jun)))
boots_mean.tb.jun
#Approximated distribution of mean
#boots_mean.tb.jun%>%
  ggplot(aes(mn)) +
  geom_density() +
  labs(title = "Approximated sampling distribution of mean res.tem.jun",
       x = "mean res.tem.jun",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.tb.jun= analysis(boots.tb.jun$splits[[12]]) #select 12 split or sample
#hist(bootstrap_sample.tb.jun$res.tem.jun)
#hist(t.blow.jun$res.tem.jun)
#For "Normal" Temperature "May"
temp.jun
t.normal.jun=temp.jun[4:5,]
mean(t.normal.jun$res.tem.jun)
library(rsample)
library(tidyverse)
samples = 300 #times
boots.tn.jun= bootstraps(t.normal.jun, times=samples)
head(boots.tn.jun)
#split into 300 times
bootstrap_sample.tn.jun= analysis(boots.tn.jun$splits[[1]]) #first split
mean(bootstrap_sample.tn.jun[,3])  #mean of first split
head(bootstrap_sample.tn.jun)
mean(analysis(boots.tn.jun$splits[[1]])$res.tem.jun) #first split mean
#all split means in the colume of mn
boots_mean.tn.jun=boots.tn.jun%>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.tem.jun)))
boots_mean.tn.jun
#Approximated distribution of mean
#boots_mean.tn.jun%>%
  ggplot(aes(mn)) +
  geom_density() +
  labs(title = "Approximated sampling distribution of mean res.tem.jun",
       x = "mean res.tem.jun",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.tn.jun= analysis(boots.tn.jun$splits[[197]]) #select 197 split or sample
#hist(bootstrap_sample.tn.jun$res.tem.jun)
#hist(t.normal.jun$res.tem.jun)
#For "Above" Temperature
temp.jun
t.above.jun=temp.jun[1:3,]
t.above.jun
mean(t.above.jun$res.tem.jun)
library(rsample)
library(tidyverse)
samples = 300 #times
boots.ta.jun= bootstraps(t.above.jun, times=samples)
head(boots.ta.jun)
#split into 300 times
bootstrap_sample.ta.jun= analysis(boots.ta.jun$splits[[1]]) #first split
mean(bootstrap_sample.ta.jun[,3])  #mean of first split
head(bootstrap_sample.ta.jun)
mean(analysis(boots.ta.jun$splits[[1]])$res.tem.jun) #first split mean
#all split means in the colume of mn
boots_mean.ta.jun=boots.ta.jun%>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.tem.jun)))
boots_mean.ta.jun
#Approximated distribution of mean
#boots_mean.ta.jun%>%
  ggplot(aes(mn)) +
  geom_density() +
  labs(title = "Approximated sampling distribution of mean res.tem.jun",
       x = "mean res.tem.jun",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.ta.jun= analysis(boots_mean.ta.jun$splits[[156]]) #select 156 split or sample
#hist(bootstrap_sample.ta.jun$res.tem.jun)
#hist(t.above.jun$res.tem.jun)
#Combine all categories of residuals into one data frame:
t.residuals.jun=rbind(bootstrap_sample.ta.jun,bootstrap_sample.tb.jun,bootstrap_sample.tn.jun)
mean(t.residuals.jun$res.tem.jun)
mean(temp.jun$res.tem.jun)
t.residuals.jun #resample residuls

#############################################
            #Precipation-3 "Jun" Ranking
#############################################
jun.pre3
#Ranking
ss=sort(jun.pre3)
ss
rank(ss)
r1=range(ss[1:2])
r1
r2=range(ss[3:4])
r2
r3=range(ss[5:7])
r3



w=function(x){
  
  df=NULL    
  for (i in 1:length(x)){
    if((x[i]>34.27)&(x[i]<64.91)){
      y=x[i]
      fac=c("Below")
      df = rbind(df,data.frame(y,fac))
    }else if((x[i]>64.90)&(x[i]<80.92)){
      y=x[i]
      fac=c("Normal")
      df = rbind(df,data.frame(y,fac))
      
    }else if((x[i]>80.91)){
      
      y=x[i]
      fac=c("Above")
      df = rbind(df,data.frame(y,fac))
      
    }
    
  }
  return(df)
}
output.pre3.jun=w(jun.pre3)
output.pre3.jun

##########################################################
#Ordering data with respect to factor Precipations-3 "Jun"
##########################################################
P.Residuals=Pfit$residuals #Residuals from Pre3 model  #Whole series

pp.residual=ts(P.Residuals,start=c(2012,1),end=c(2018,12),frequency=12)

res.pre.jun=subset(pp.residual,cycle(pp.residual)==6)
res.pre.jun
year=rep(c(2012:2018), each = 1)

app=data.frame(year,output.pre3.jun,res.pre.jun)
df.pre3.jun<- app[order(app$fac),]
df.pre3.jun #Order wise data

##########################################################
#Conditional Residuls resampling of Precipation-3 "Jun"
##########################################################
df.pre3.jun #Order wise data

pre.3m=df.pre3.jun [,-2]
pre.3m
#For "Below" Precipation3 
p.blow.jun =pre.3m[6:7,]
mean(p.blow.jun $res.pre.jun )
library(rsample)
library(tidyverse)
samples = 500 #times
boots.pb.jun = bootstraps(p.blow.jun , times=samples)
head(boots.pb.jun )
#split into 300 times
bootstrap_sample.pb.jun = analysis(boots.pb.jun$splits[[1]]) #first split
mean(bootstrap_sample.pb.jun[,3])  #mean of first split
head(bootstrap_sample.pb.jun)
mean(analysis(boots.pb.jun$splits[[1]])$res.pre.jun) #first split mean
#all split means in the colume of mn
boots_mean.pb.jun = boots.pb.jun%>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.pre.jun)))
boots_mean.pb.jun
#Approximated distribution of mean
#boots_mean.pb.jun%>%
  ggplot(aes(mn)) +
  geom_density() +
  labs(title = "Approximated sampling distribution of mean res.pre.jun",
       x = "mean res.pre.jun",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.pb.jun = analysis(boots.pb.jun$splits[[27]]) #select 27 split or sample
#hist(bootstrap_sample.pb.jun$res.pre.jun)
#hist(p.blow.jun$res.pre.jun)

#For "Normal" Preciaption3 
pre.3m
p.normal.jun=pre.3m[4:5,]
mean(p.normal.jun$res.pre.jun)
library(rsample)
library(tidyverse)
samples = 500 #times
boots.pn.jun = bootstraps(p.normal.jun, times=samples)
head(boots.pn.jun)
#split into 300 times
bootstrap_sample.pn.jun = analysis(boots.pn.jun$splits[[1]]) #first split
mean(bootstrap_sample.pn.jun[,3])  #mean of first split
head(bootstrap_sample.pn.jun)
mean(analysis(boots.pn.jun$splits[[1]])$res.pre.jun) #first split mean
#all split means in the colume of mn
boots_mean.pn.jun= boots.pn.jun%>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.pre.jun)))
boots_mean.pn.jun
#Approximated distribution of mean
#boots_mean.pn.jun%>%
  ggplot(aes(mn)) +
  geom_density() +
  labs(title = "Approximated sampling distribution of mean res.pre.jun",
       x = "mean res.pre.jun",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.pn.jun= analysis(boots.pn.jun$splits[[29]]) #select 29 split or sample
#hist(bootstrap_sample.pn.jun$res.pre.jun)
#hist(p.normal.jun$res.pre.jun)

#For "Above" Preciaption3 
pre.3m
p.above.jun=pre.3m[1:3,]
mean(p.normal.jun$res.pre.jun)
library(rsample)
library(tidyverse)
samples = 500 #times
boots.pa.jun= bootstraps(p.above.jun, times=samples)
head(boots.pa.jun)
#split into 500 times
bootstrap_sample.pa.jun= analysis(boots.pa.jun$splits[[1]]) #first split
mean(bootstrap_sample.pa.jun[,3])  #mean of first split
head(bootstrap_sample.pa.jun)
mean(analysis(boots.pa.jun$splits[[1]])$res.pre.jun) #first split mean
#all split means in the colume of mn
boots_mean.pa.jun= boots.pa.jun%>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.pre.jun)))
boots_mean.pa.jun
#Approximated distribution of mean
#boots_mean.pa.jun%>%
  ggplot(aes(mn)) +
  geom_density() +
  labs(title = "Approximated sampling distribution of mean res.pre.jun",
       x = "mean res.pre.jun",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.pa.jun= analysis(boots.pa.jun$splits[[77]]) #select 77 split or sample
#hist(bootstrap_sample.pa.jun$res.pre.jun)
#hist(p.above.jun$res.pre.jun)
#Combine all categories of residuals into one data frame:
p.residuals.jun=rbind(bootstrap_sample.pn.jun,bootstrap_sample.pa.jun,bootstrap_sample.pb.jun)
mean(p.residuals.jun$res.pre.jun)
mean(pre.3m$res.pre.jun)
p.residuals.jun #resample residuls
#############################################################
         #End for Jun
#############################################################

#############################################################
         #Work for July
#############################################################
#Rank the observed of target month "MAY" temperature and pre3 into three categories 
jul.tem

length(jul.tem)
ss=sort(jul.tem)
ss
rank(ss)
r1=range(ss[1:2])
r1
r2=range(ss[3:4])
r2
r3=range(ss[5:6])
r3
w=function(x){
  
  df=NULL    
  for (i in 1:length(x)){
    if((x[i]>14.41)&(x[i]<18.11)){
      y=x[i]
      fac=c("Below")
      df = rbind(df,data.frame(y,fac))
    }else if((x[i]>18.10)&(x[i]<22.06)){
      y=x[i]
      fac=c("Normal")
      df = rbind(df,data.frame(y,fac))
      
    }else if((x[i]>22.05)){
      
      y=x[i]
      fac=c("Above")
      df = rbind(df,data.frame(y,fac))
      
    }
    
  }
  return(df)
}

output.tem.jul=w(jul.tem)
output.tem.jul
#########################################################
#Ordering data with respect to factor TEMPERATURE "July"
#########################################################
T.Residuals=Tfit$residuals #Residuals from temperature model  #Whole series
length(T.Residuals)
tt.residual=ts(T.Residuals,start=c(2012,1),end=c(2018,12),frequency=12)
#tt.residual=ts(T.Residuals,start=c(1990,1),end=c(2015,12),frequency=12)
res.tem.jul=subset(tt.residual,cycle(tt.residual)==7)
res.tem.jul
year=rep(c(2012:2018), each = 1)
app=data.frame(year,output.tem.jul,res.tem.jul)
df.tem.jul<- app[order(app$fac),]
df.tem.jul#Order wise data


##########################################################
#Conditional Residuls resampling of temperature "July"
##########################################################
df.tem.jul
temp.jul=df.tem.jul[,-2]
#For "Below Temperature" 
temp.jul
t.blow.jul=temp.jul[6:7,]
mean(t.blow.jul$res.tem.jul)
library(rsample)
library(tidyverse)
samples = 300 #times
boots.tb.jul= bootstraps(t.blow.jul, times=samples)
head(boots.tb.jul)
#split into 300 times
bootstrap_sample.tb.jul= analysis(boots.tb.jul$splits[[1]]) #first split
mean(bootstrap_sample.tb.jul[,3])  #mean of first split
head(bootstrap_sample.tb.jul)
mean(analysis(boots.tb.jul$splits[[1]])$res.tem.jul) #first split mean
#all split means in the colume of mn
boots_mean.tb.jul= boots.tb.jul%>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.tem.jul)))
boots_mean.tb.jul
#Approximated distribution of mean
#boots_mean.tb.jul%>%
  ggplot(aes(mn)) +
  geom_density() +
  labs(title = "Approximated sampling distribution of mean res.tem.jul",
       x = "mean res.tem.jul",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.tb.jul= analysis(boots.tb.jul$splits[[12]]) #select 12 split or sample
#hist(bootstrap_sample.tb.jul$res.tem.jul)
#hist(t.blow.jun$res.tem.jul)

#For "Normal" Temperature "May"
temp.jul
t.normal.jul=temp.jul[4:5,]
mean(t.normal.jul$res.tem.jul)
library(rsample)
library(tidyverse)
samples = 300 #times
boots.tn.jul= bootstraps(t.normal.jul, times=samples)
head(boots.tn.jul)
#split into 300 times
bootstrap_sample.tn.jul= analysis(boots.tn.jul$splits[[1]]) #first split
mean(bootstrap_sample.tn.jul[,3])  #mean of first split
head(bootstrap_sample.tn.jul)
mean(analysis(boots.tn.jul$splits[[1]])$res.tem.jul) #first split mean
#all split means in the colume of mn
boots_mean.tn.jul=boots.tn.jul %>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.tem.jul)))
boots_mean.tn.jul
#Approximated distribution of mean
#boots_mean.tn.jul%>%
  ggplot(aes(mn)) +
  geom_density() +
  labs(title = "Approximated sampling distribution of mean res.tem.jul",
       x = "mean res.tem.jul",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.tn.jul= analysis(boots.tn.jul$splits[[197]]) #select 197 split or sample
#hist(bootstrap_sample.tn.jul$res.tem.jul)
#hist(t.normal.jul$res.tem.jul)
#For "Above" Temperature
temp.jul
t.above.jul=temp.jul[1:3,]
t.above.jul
mean(t.above.jul$res.tem.jul)
library(rsample)
library(tidyverse)
samples = 300 #times
boots.ta.jul= bootstraps(t.above.jul, times=samples)
head(boots.ta.jul)
#split into 300 times
bootstrap_sample.ta.jul= analysis(boots.ta.jul$splits[[1]]) #first split
mean(bootstrap_sample.ta.jul[,3])  #mean of first split
head(bootstrap_sample.ta.jul)
mean(analysis(boots.ta.jul$splits[[1]])$res.tem.jul) #first split mean
#all split means in the colume of mn
boots_mean.ta.jul=boots.ta.jul%>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.tem.jul)))
boots_mean.ta.jul
#Approximated distribution of mean
#boots_mean.ta.jul%>%
  ggplot(aes(mn)) +
  geom_density() +
  labs(title = "Approximated sampling distribution of mean res.tem.jul",
       x = "mean res.tem.jul",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.ta.jul= analysis(boots_mean.ta.jul$splits[[156]]) #select 156 split or sample
#hist(bootstrap_sample.ta.jul$res.tem.jul)
#hist(t.above.jul$res.tem.jul)
#Combine all categories of residuals into one data frame:
t.residuals.jul=rbind(bootstrap_sample.ta.jul,bootstrap_sample.tb.jul,bootstrap_sample.tn.jul)
mean(t.residuals.jul$res.tem.jul)
mean(temp.jul$res.tem.jul)
t.residuals.jul #resample residuls


#############################################
            #Precipation-3 "July" Ranking
#############################################
jul.pre3
#Ranking
ss=sort(jul.pre3)
ss
rank(ss)
r1=range(ss[1:2])
r1
r2=range(ss[3:4])
r2
r3=range(ss[5:7])
r3

w=function(x){
  
  df=NULL    
  for (i in 1:length(x)){
    if((x[i]>25.94)&(x[i]<47.15)){
      y=x[i]
      fac=c("Below")
      df = rbind(df,data.frame(y,fac))
    }else if((x[i]>47.14)&(x[i]<92.3)){
      y=x[i]
      fac=c("Normal")
      df = rbind(df,data.frame(y,fac))
      
    }else if((x[i]>92.2)){
       y=x[i]
      fac=c("Above")
      df = rbind(df,data.frame(y,fac))
      
    }
    
  }
  return(df)
}

output.pre3.jul=w(jul.pre3)
output.pre3.jul

#########################################################
#Ordering data with respect to factor Precipations-3 "July"
#########################################################
P.Residuals=Pfit$residuals #Residuals from Pre3 model  #Whole series

pp.residual=ts(P.Residuals,start=c(2012,1),end=c(2018,12),frequency=12)
res.pre.jul=subset(pp.residual,cycle(pp.residual)==7)
res.pre.jul
year=rep(c(2012:2018), each = 1)
app=data.frame(year,output.pre3.jul,res.pre.jul)
df.pre3.jul<- app[order(app$fac),]
df.pre3.jul#Order wise data

##########################################################
#Conditional Residuls resampling of Precipation-3 "July"
##########################################################
df.pre3.jul #Order wise data

pre.3j=df.pre3.jul [,-2]
pre.3j
#For "Below" Precipation3 
p.blow.jul =pre.3j[6:7,]
mean(p.blow.jul $res.pre.jul )
library(rsample)
library(tidyverse)
samples = 500 #times
boots.pb.jul = bootstraps(p.blow.jul , times=samples)
head(boots.pb.jul )
#split into 300 times
bootstrap_sample.pb.jul = analysis(boots.pb.jul $splits[[1]]) #first split
mean(bootstrap_sample.pb.jul [,3])  #mean of first split
head(bootstrap_sample.pb.jul )
mean(analysis(boots.pb.jul $splits[[1]])$res.pre.jul ) #first split mean
#all split means in the colume of mn
boots_mean.pb.jul = boots.pb.jul %>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.pre.jul )))
boots_mean.pb.jul 
#Approximated distribution of mean
#boots_mean.pb.jul %>%
  ggplot(aes(mn)) +
  geom_density() +
  labs(title = "Approximated sampling distribution of mean res.pre.jul ",
       x = "mean res.pre.jul ",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.pb.jul = analysis(boots.pb.jul $splits[[27]]) #select 27 split or sample
#hist(bootstrap_sample.pb.jul $res.pre.jul )
#hist(p.blow.jul $res.pre.jul)

#For "Normal" Preciaption3 
pre.3j
p.normal.jul =pre.3j[4:5,]
mean(p.normal.jul $res.pre.jul)
library(rsample)
library(tidyverse)
samples = 500 #times
boots.pn.jul= bootstraps(p.normal.jul, times=samples)
head(boots.pn.jul)
#split into 300 times
bootstrap_sample.pn.jul= analysis(boots.pn.jul$splits[[1]]) #first split
mean(bootstrap_sample.pn.jul[,3])  #mean of first split
head(bootstrap_sample.pn.jul)
mean(analysis(boots.pn.jul$splits[[1]])$res.pre.jul) #first split mean
#all split means in the colume of mn
boots_mean.pn.jul= boots.pn.jul%>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.pre.jul)))
boots_mean.pn.jul
#Approximated distribution of mean
#boots_mean.pn.jul%>%
  ggplot(aes(mn)) +
  geom_density() +
  labs(title = "Approximated sampling distribution of mean res.pre.jul",
       x = "mean res.pre.jul",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.pn.jul= analysis(boots.pn.jul$splits[[29]]) #select 29 split or sample
#hist(bootstrap_sample.pn.jul$res.pre.jul)
#hist(p.normal.jul$res.pre.jul)

#For "Above" Preciaption3 
pre.3j
p.above.jul=pre.3j[1:3,]
mean(p.normal.jul$res.pre.jul)
library(rsample)
library(tidyverse)
samples = 500 #times
boots.pa.jul= bootstraps(p.above.jul, times=samples)
head(boots.pa.jul)
#split into 500 times
bootstrap_sample.pa.jul= analysis(boots.pa.jul$splits[[1]]) #first split
mean(bootstrap_sample.pa.jul[,3])  #mean of first split
head(bootstrap_sample.pa.jul)
mean(analysis(boots.pa.jul$splits[[1]])$res.pre.jul) #first split mean
#all split means in the colume of mn
boots_mean.pa.jul= boots.pa.jul%>%
  mutate(mn = map_dbl(splits, ~ mean(analysis(.x)$res.pre.jul)))
boots_mean.pa.jul
#Approximated distribution of mean
#boots_mean.pa.jul%>%
  ggplot(aes(mn)) +
  geom_density() +
  labs(title = "Approximated sampling distribution of mean res.pre.jul",
       x = "mean res.pre.jul",
       y = "frequency") +
  theme_minimal()
#Select only one split 
#Randomly select split
bootstrap_sample.pa.jul= analysis(boots.pa.jul$splits[[77]]) #select 77 split or sample
#hist(bootstrap_sample.pa.jul$res.pre.jul)
#hist(p.above.jul$res.pre.jul)
#Combine all categories of residuals into one data frame:
p.residuals.jul=rbind(bootstrap_sample.pn.jul,bootstrap_sample.pa.jul,bootstrap_sample.pb.jul)
mean(p.residuals.jul$res.pre.jul)
mean(pre.3j$res.pre.jul)
p.residuals.jul #resample residuls pre3 of jun
########################################################
               #End for July
########################################################


####################################################################
#Observe data frame test data set 2010-2015 for comparison perpose
####################################################################
#Test data 2012-2018 From SPI3
Test.SP3=spi3[265:348] #testing set
#COnvert into time series data and select data:


Observe.Months=ts(Test.SP3,start=c(2012,1),end=c(2018,12),frequency=12)
#Get all data from class of "ts"

#Target Months
Obs.May=subset(Observe.Months, cycle(Observe.Months) == 5)
Obs.May
Obs.Jun=subset(Observe.Months, cycle(Observe.Months) == 6)
Obs.Jun
Obs.Jul=subset(Observe.Months, cycle(Observe.Months) == 7)
Obs.Jul

####################################################################

########################################################
            #Ensemble Forecast of Months:
########################################################

#Ensemble Forecast of "May"
#May.ensemble=y00.may+ee1.may+t.residuals.may$res.tem.may+p.residuals.may$res.pre.may
#May.ensemble
#range(scale(May.ensemble))
#range(May.ensemble)
#may.for=scale(May.ensemble)
#may.for=as.numeric(may.for)  #Ensemble Forecast
#range(may.for)

May.Int.RS=scale(y00.may)
May.Res.RS=scale(ee1.may+t.residuals.may$res.tem.may+p.residuals.may$res.pre.may)
May.Res.RS=as.numeric(May.Res.RS)


#Ensemble Forecast of "Jun"
#Jun.ensemble=y00.jun+ee1.jun+t.residuals.jun$res.tem.jun+p.residuals.jun$res.pre.jun
#Jun.ensemble
#range(scale(Jun.ensemble))
#range(Jun.ensemble)
#jun.for=scale(Jun.ensemble)
#jun.for=as.numeric(jun.for)
#range(jun.for)

Jun.Int.RS=scale(y00.jun)
Jun.Res.RS=scale(ee1.jun+t.residuals.jun$res.tem.jun+p.residuals.jun$res.pre.jun)
Jun.Res.RS=as.numeric(Jun.Res.RS)


#Ensemble Forecast of "Jul"
#Jul.ensemble=y00.jul+ee1.jun+t.residuals.jul$res.tem.jul+p.residuals.jul$res.pre.jul
#Jul.ensemble
#range(scale(Jul.ensemble))
#range(Jul.ensemble)
#jul.for=scale(Jul.ensemble)
#jul.for=as.numeric(jul.for)

Jul.Int.RS=scale(y00.jul)
Jul.Res.RS=scale(ee1.jul+t.residuals.jul$res.tem.jul+p.residuals.jul$res.pre.jul)
Jul.Res.RS=as.numeric(Jul.Res.RS)

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

Modelperformance(may.for,Obs.May)
Modelperformance(jun.for,Obs.Jun)
Modelperformance(jul.for,Obs.Jul)


###################################################################
#Equal Weighting Method:
#make ensmeble member on the basis of years in the series:

###################
#Equal Ensmebling:
###################

######
#May
######
set.seed(333)
Resamples <- lapply(1:21, function(i) sample(May.Int.RS  , replace=FALSE))
Resamples
EE1=Resamples[1]
EE2=Resamples[2]
EE3=Resamples[3]
EE4=Resamples[4]
EE5=Resamples[5]
EE6=Resamples[6]
EE7=Resamples[7]

Ensmeble.May=cbind.data.frame(EE1,EE2,EE3,EE4,EE5,EE6,EE7)
colnames(Ensmeble.May)=c("member1","member2","member3","member4","member5","member6","member7")  
Ensmeble.May
N=7
W=1/N
W

May.EEF=Ensmeble.May+May.Res.RS+W
May.EEDP<-rowMeans(May.EEF[,1:7]) 
range(May.EEDP)  
Modelperformance(May.EEDP,Obs.May)
#boxplot(May.EEF)

#####
#Jun
#####
set.seed(333)
Resamples <- lapply(1:21, function(i) sample(Jun.Int.RS  , replace=FALSE))
Resamples
EE1=Resamples[1]
EE2=Resamples[2]
EE3=Resamples[3]
EE4=Resamples[4]
EE5=Resamples[5]
EE6=Resamples[6]
EE7=Resamples[7]

Ensmeble.Jun=cbind.data.frame(EE1,EE2,EE3,EE4,EE5,EE6,EE7)
colnames(Ensmeble.Jun)=c("member1","member2","member3","member4","member5","member6","member7")  
Ensmeble.Jun
N=7
W=1/N
W


Jun.EEF=Ensmeble.Jun+Jun.Res.RS+W
Jun.EEDP<-rowMeans(Jun.EEF[,1:7]) 
range(Jun.EEDP)  
Modelperformance(Jun.EEDP,Obs.Jun)
#boxplot(Jun.EEF)
#####
#Jul
######
set.seed(333)
Resamples <- lapply(1:21, function(i) sample(Jul.Int.RS  , replace=FALSE))
Resamples
EE1=Resamples[1]
EE2=Resamples[2]
EE3=Resamples[3]
EE4=Resamples[4]
EE5=Resamples[5]
EE6=Resamples[6]
EE7=Resamples[7]

Ensmeble.Jul=cbind.data.frame(EE1,EE2,EE3,EE4,EE5,EE6,EE7)
colnames(Ensmeble.Jul)=c("member1","member2","member3","member4","member5","member6""member7")  
Ensmeble.Jul
N=7
W=1/N
W

Jul.EEF=Ensmeble.Jul+Jul.Res.RS+W
Jul.EEDP<-rowMeans(Jul.EEF[,1:7]) 
range(Jul.EEDP)  
Modelperformance(Jul.EEDP,Obs.Jul)
#boxplot(Jun.EEF)

#EEDP Burkhan
Modelperformance(May.EEDP,Obs.May)
Modelperformance(Jun.EEDP,Obs.Jun)
Modelperformance(Jul.EEDP,Obs.Jul) 
##############################################################################################################
#############################
#Cliate Predictor 
###############################
#Climate predictor for obs data
#fix(data1)
may.pre.obs=May.PRE.obs
jun.pre.obs=Jun.PRE.obs
jul.pre.obs=Jul.PRE.obs


May.mix.obs=data.frame(may.pre.obs,May.AM.obs,May.NTA.obs,May.TNAI.obs)
May.Model.obs=lm(may.pre.obs~May.AM.obs+May.NTA.obs+May.TNAI.obs, data=May.mix.obs)
summary(May.Model.obs)
Integrated.Climate.May.Obs=May.Model.obs$fitted
Integrated.Climate.May.Obs=as.numeric(Integrated.Climate.May.Obs)

#claimate predictor for predict data

May.mix.pred=data.frame(May.PRE.pred,May.AM.pred,May.NTA.pred,May.TNAI.pred)
May.Model.pred=lm(May.PRE.pred~May.AM.pred+May.NTA.pred+May.TNAI.pred, data=May.mix.pred)
summary(May.Model.pred)
Integrated.Climate.May.Pred=May.Model.pred$fitted
Integrated.Climate.May.Pred=as.numeric(Integrated.Climate.May.Pred)
Integrated.Climate.forcast.May=Integrated.Climate.May.Pred
################################################################################################################

#Weighted Ensemble Model

#Weighted Ensemble Model


set.seed(300)
C.sample <- lapply(1:21, function(i) sample(Integrated.Climate.May.Obs, replace=FALSE))
C.sample
cc1=C.sample[1]
cc2=C.sample[2]
cc3=C.sample[3]
cc4=C.sample[4]
cc5=C.sample[5]
cc6=C.sample[6]
cc7=C.sample[7]

Climate.Rsample=cbind.data.frame(cc1,cc2,cc3,cc4,cc5,cc6,cc7)
colnames(Climate.Rsample)=c("theta1","theta2","theta3","theta4","theta5","theta6","theta7")
Climate.Rsample
f1=abs(Climate.Rsample$theta1-Integrated.Climate.forcast.May)
f2=abs(Climate.Rsample$theta2-Integrated.Climate.forcast.May)
f3=abs(Climate.Rsample$theta3-Integrated.Climate.forcast.May)
f4=abs(Climate.Rsample$theta4-Integrated.Climate.forcast.May)
f5=abs(Climate.Rsample$theta5-Integrated.Climate.forcast.May)
f6=abs(Climate.Rsample$theta6-Integrated.Climate.forcast.May)
f7=abs(Climate.Rsample$theta7-Integrated.Climate.forcast.May)

#For Kernal Function:
#install.packages("extremefit")
library("extremefit")
#Truncated gussain function

Sum.C=(abs(f1)+abs(f2)+abs(f3)+abs(f4)+abs(f5)+abs(f6)+abs(f7))


sigma1=sd(f1)
Kx1=TruncGauss.kernel(f1, sigma1)
Kx1
WW1=Kx1*(abs(f1))/(Kx1*Sum.C)
WW1[which(!is.finite(WW1))] <- 0.000000 #Replace Inf
WW1=data.frame(WW1)
WW1=rapply(WW1, f=function(x) ifelse(is.nan(x),0,x), how="replace" ) #Replace NaN
WW1=as.matrix(WW1)
Weight1=as.numeric(WW1)



sigma2=sd(f2)
Kx2=TruncGauss.kernel(f2, sigma2)
Kx2
WW2=Kx2*(abs(f2))/(Kx2*Sum.C)
WW2[which(!is.finite(WW2))] <- 0.000000 #Replace Inf
WW2=data.frame(WW2)
WW2=rapply(WW2, f=function(x) ifelse(is.nan(x),0,x), how="replace" ) #Replace NaN
WW2=as.matrix(WW2)
Weight2=as.numeric(WW2)

sigma3=sd(f3)
Kx3=TruncGauss.kernel(f3, sigma3)
Kx3
WW3=Kx3*(abs(f3))/(Kx3*Sum.C)
WW3[which(!is.finite(WW3))] <- 0.000000 #Replace Inf
WW3=data.frame(WW3)
WW3=rapply(WW3, f=function(x) ifelse(is.nan(x),0,x), how="replace" ) #Replace NaN
WW3=as.matrix(WW3)
Weight3=as.numeric(WW3)

sigma4=sd(f4)
Kx4=TruncGauss.kernel(f4, sigma4)
Kx4
WW4=Kx4*(abs(f4))/(Kx4*Sum.C)
WW4[which(!is.finite(WW4))] <- 0.000000 #Replace Inf
WW4=data.frame(WW4)
WW4=rapply(WW4, f=function(x) ifelse(is.nan(x),0,x), how="replace" ) #Replace NaN
WW4=as.matrix(WW4)
Weight4=as.numeric(WW4)

sigma5=sd(f5)
Kx5=TruncGauss.kernel(f5, sigma5)
Kx5
WW5=Kx5*(abs(f5))/(Kx5*Sum.C)
WW5[which(!is.finite(WW5))] <- 0.000000 #Replace Inf
WW5=data.frame(WW5)
WW5=rapply(WW5, f=function(x) ifelse(is.nan(x),0,x), how="replace" ) #Replace NaN
WW5=as.matrix(WW5)
Weight5=as.numeric(WW5)

sigma6=sd(f6)
Kx6=TruncGauss.kernel(f6, sigma6)
Kx6
WW6=Kx6*(abs(f6))/(Kx6*Sum.C)
WW6[which(!is.finite(WW6))] <- 0.000000 #Replace Inf
WW6=data.frame(WW6)
WW6=rapply(WW6, f=function(x) ifelse(is.nan(x),0,x), how="replace" ) #Replace NaN
WW6=as.matrix(WW6)
Weight6=as.numeric(WW6)

sigma7=sd(f7)
Kx7=TruncGauss.kernel(f7, sigma7)
Kx7
WW7=Kx6*(abs(f7))/(Kx7*Sum.C)
WW7[which(!is.finite(WW7))] <- 0.000000 #Replace Inf
WW7=data.frame(WW7)
WW7=rapply(WW7, f=function(x) ifelse(is.nan(x),0,x), how="replace" ) #Replace NaN
WW7=as.matrix(WW7)
Weight7=as.numeric(WW7)

#Give Weight Each Ensemble Memebers

Ens1=Ensmeble.May$member1+Weight1
Ens2=Ensmeble.May$member2+Weight2
Ens3=Ensmeble.May$member3+Weight3
Ens4=Ensmeble.May$member4+Weight4
Ens5=Ensmeble.May$member5+Weight5
Ens6=Ensmeble.May$member6+Weight6
Ens7=Ensmeble.May$member7+Weight7
W.Enemble.May=cbind.data.frame(Ens1,Ens2,Ens3,Ens4,Ens5,Ens6,Ens7)

May.WEF=Ensmeble.May+May.Res.RS
May.WEDP<-rowMeans(May.WEF[,1:7]) 
range(May.WEDP)  
Modelperformance(May.WEDP,Obs.May)
########################################################################################

Weighted bootstrap procedure:

ww1=BT.Int.May$member1/BT.ICI.May$theta1
ww2=BT.Int.May$member2/BT.ICI.May$theta2
ww3=BT.Int.May$member3/BT.ICI.May$theta3
ww4=BT.Int.May$member4/BT.ICI.May$theta4
ww5=BT.Int.May$member5/BT.ICI.May$theta5
ww6=BT.Int.May$member6/BT.ICI.May$theta6
ww7=BT.Int.May$member7/BT.ICI.May$theta7

QQ=(ww1+ww2+ww3+ww4+ww5+ww6+ww7)
q1=ww1/QQ
q2=ww2/QQ
q3=ww3/QQ
q4=ww4/QQ
q5=ww5/QQ
q6=ww6/QQ
q7=ww7/QQ

#Give Bayesion Weight Each Ensemble Memebers
May.Res.RS=as.numeric(May.Res.RS)

Ensb1=BT.Int.May$member1+q1+May.Res.RS
Ensb2=BT.Int.May$member2+q2+May.Res.RS
Ensb3=BT.Int.May$member3+q3+May.Res.RS
Ensb4=BT.Int.May$member4+q4+May.Res.RS
Ensb5=BT.Int.May$member5+q5+May.Res.RS
Ensb6=BT.Int.May$member6+q6+May.Res.RS
Ensb7=BT.Int.May$member7+q7+May.Res.RS

bayes.Enemble.May=cbind.data.frame(Ensb1,Ensb2,Ensb3,Ensb4,Ensb5,Ensb6,Ensb7)

May.Ensemble.RMB<-rowMeans(bayes.Enemble.May[,1:7])                                                
May.Ensemble.RMB
Modelperformance(May.Ensemble.RMB,Obs.May)
########################################################################################
Conditional Weighting Ensemble Model

#load may precipiation data obs
PRE.data=read.csv(file.choose())
May.PRE=PRE.data$May[23:29]
Jun.PRE=PRE.data$Jun[23:29]
Jul.PRE=PRE.data$May[23:29]
#########
#FOR May
#########

#By Using Couple Function:

#install.packages("copula")
library("copula")
#Fit log normal distribution on both variable
library(MASS)
#April Precipation
LN.Pre <- fitdistr(May.PRE, "lognormal") #parameter Estimate
LN.Pre
Fit.LN.May.P=dlnorm(May.PRE, meanlog =2.2407092 , sdlog =1.4177163,log = FALSE)
CFit.LN.May.P=plnorm(May.PRE,  meanlog =2.2407092 , sdlog =1.4177163, log = FALSE)

#Climate Index fit Ditribution
LN.CI.May <- fitdistr(Integrated.Climate.May.Obs, "lognormal") #parameter Estimate
LN.CI.May
Fit.LN.May.CI=dlnorm(Integrated.Climate.May.Obs, meanlog =3.4541949 , sdlog =1.2763469, log = FALSE)
Fit.LN.May.CI
CFit.LN.May.CI=plnorm(Integrated.Climate.May.Obs, meanlog =3.4541949 , sdlog =1.2763469, log = FALSE)
CFit.LN.May.CI

#Fit Different Family of Coupla Functions:

#check approprate coupla
var_a <- pobs(CFit.LN.May.P)
var_b <- pobs(CFit.LN.May.CI)
cor(var_a,var_b)
library(VineCopula)
#automatic select family on the basis of bic value
Apro.family=BiCopSelect(var_a, var_b, familyset =6, selectioncrit = "BIC",
  indeptest = FALSE, level = 0.05, weights = NA, rotations = TRUE,
  se = FALSE, presel = TRUE, method = "mle")

summary(Apro.family) #frank is aaproprite

#Kendall  log  #aic  #bic 

#BIC
#frank 5 
#gumbel 4 
#clayton 3 
#joe 

#Using Joe copula family

joe.model <- joeCopula(param =1.23, dim = 2)
mydata.may=data.frame(CFit.LN.May.P,CFit.LN.May.CI)
m.may <- pobs(as.matrix(mydata.may))

joe.fit <- fitCopula(joe.model,m.may , method = 'ml')
coef(joe.fit) #alpha 1.230574
str(joe.fit)
# Check Kendall's tau value for the joe copula
tau(joeCopula(param =1.23)) #tau 0.1159444

#Building the bivariate distribution using the clayton Family

margins = c("lnorm","lnorm")
paramMargins = list(list(meanlog =2.2407092 , sdlog =1.4177163), list(meanlog =3.4541949 , sdlog =1.2763469))

Joe.dist <- mvdc(joeCopula(param =1.23, dim = 2),
 margins, paramMargins,marginsIdentical =FALSE)
Joe.dist

Generate random sample observations from the multivariate distribution
vc <- rMvdc(7, Joe.dist)
# Compute the density
Cpdf_mvd <- dMvdc(vc, Joe.dist)
# Compute the CDF
Ccdf_mvd <- pMvdc(vc,Joe.dist)

Joe.CD=Cpdf_mvd /(Fit.LN.May.CI)
Joe.CD
hist(Joe.CD)
#Check Distribution of conditional density
library("fitdistrplus")
library("propagate")
joe.f=fitDistr(Joe.CD, nbin = 100, weights = FALSE,verbose = TRUE,  brute = c("fast", "slow"), 
            plot = c("hist", "qq"),  distsel = NULL)
joe.f #johnson sb #bic 536.43
#check through easy fit software
Joe.CD.Easy=as.matrix(Joe.CD) #
Joe.CD.Easy 

#install.packages("PearsonDS")
library(PearsonDS)

MSC  <- pearsonMSC(Gumbel.CD,control=list(iter.max=1e5,eval.max=1e5))
## log-Likelihood values for all distribution sub-classes
print(MSC$logLik)
## Values for all MSCs and distribution sub-classes
print(MSC$MSCs)  #BIC -1.4530094  AIC -0.8282878

#broom::glance(MASS::fitdistr(Frank.CD,"weibull"))  

library("stats")
#log normal is approporate 
#BIC 
#May.Con.Fit=dlnorm(JOE.CD, meanlog =-1.634 , sdlog =2.8735, log = FALSE)
#CMay.Con.Fit=plnorm(JOE.CD, meanlog =-1.634 , sdlog =2.8735, lower.tail = TRUE, log.p = FALSE)
plot(JOE.CD,CMay.Con.Fit)


#GEV is apprrprate
library(fExtremes)
#May.Con.Fit=dgev(Frank.CD, xi =0.61254, mu =0.05452, beta =0.09311, log = FALSE)
#CMay.Con.Fit=pgev(Frank.CD, xi =0.61254, mu =0.05452, beta =0.09311, lower.tail = TRUE)


#Frechet Distribution
library(evd)
#May.Con.Fit=dfrechet(Joe.CD, loc=0, scale=0.00218, shape=0.33393, log = FALSE)
#CMay.Con.Fit=pfrechet(Joe.CD, loc=0, scale=0.00218, shape=0.33393, lower.tail = TRUE) 
plot(Joe.CD,CMay.Con.Fit)

#Gamma Distribution 
library(stats)
May.Con.Fit=dgamma(Joe.CD, shape=0.65139, scale =0.08513, log = FALSE)
CMay.Con.Fit=pgamma(Joe.CD, shape=0.65139, scale =0.08513, lower.tail = TRUE,
       log.p = FALSE)

m=0
length(CMay.Con.Fit)
n=length(CMay.Con.Fit)
q=m/n  # m is total number of zero present in present time series
#n is total number of observations # q is probablity of zero.
q
Gx=CMay.Con.Fit
Gx
Hx=q+(1-q)*Gx
Hx
range(Hx)
#SPI Loop
sp=function(Hx){
  list=c()
  c0 = 2.515517
  c1 = 0.802583
  c2 = 0.010328
  d1 = 1.432788
  d2 = 0.189269
  d3 = 0.001308 
  for (i in 1:length(Hx)) {
    
    if((Hx[i]>0)&(Hx[i]<0.5)){
      
      t1=sqrt(log(((1/Hx[i])^2)))
      m1=((c0+c1*(t1)+c2*(t1^2))/(1+d1*(t1)+d2*(t1^2)+d3*(t1^3)))
      spiz=-(t1-m1)
      list[[i]]=spiz
      
    } else if((Hx[i]>0.5)&(Hx[i]<1.0)){
      
      t2=sqrt(log((1/(1-Hx[i])^2)))
      m2=((c0+c1*(t2)+c2*(t2^2))/(1+d1*(t2)+d2*(t2^2)+d3*(t2^3)))
      spiz=+(t2-m2)
      list[[i]]=spiz
    }
    
    
  } 
  
  print(list)
}

May.Con.Pre=sp(Hx)
May.Con.Pre
Modelperformance(May.Con.Pre,Obs.May)

set.seed(1000)
resamples =lapply(1:1000,function(i) sample(May.Con.Pre, replace=FALSE))
Conditional.RS.May<- data.frame(matrix(unlist(resamples), nrow=7, byrow=F),stringsAsFactors=FALSE)
Conditional.RS.May

May.CEDP.RM<-rowMeans(Conditional.RS.May[,1:1000])
Modelperformance(May.CEDP.RM,Obs.May)