# R language code for copula # Data import dta<-read.csv("E:\\research\\drought\\copula\\SRI\\for R\\drought events.csv") head(dta) D<-dta$duration S<-dta$severity I<-dta$intensity library(copula) # dependence estimate for drought variables cor(D,S,method="pearson") cor(D,I,method="pearson") cor(S,I,method="pearson") cor(D,S,method="kendall") cor(D,I,method="kendall") cor(S,I,method="kendall") # dependence estimate for transformed uniform drought variables Ud<-pweibull(D,shape=1.9456,scale=3.616) Us<-pgamma(S,shape=2.5196,rate=0.6482) Ui<-pgamma(I,shape=113.6474,rate=97.206) cor(Ud,Us,method="pearson") cor(Ud,Ui,method="pearson") cor(Us,Ui,method="pearson") cor(Ud,Us,method="kendall") cor(Ud,Ui,method="kendall") cor(Us,Ui,method="kendall") # find the parameters (0.916,0.679,0.78) estimated by Kendall's tau ############ copula contour plot ################ #----------------Ud,Us------------------ myCop.t1<-tCopula(dim=2,param=0.9955,df=1.4468e7) myCop.g1<-gumbelCopula(13.724) myMvd1<-mvdc(copula=myCop.t1,margins=c("weibull","gamma"), paramMargins=list(list(shape=1.9456,scale=3.616), list(shape=2.5196,rate=0.6482))) myMvd2<-mvdc(copula=myCop.g1,margins=c("weibull","gamma"), paramMargins=list(list(shape=1.9456,scale=3.616), list(shape=2.5196,rate=0.6482))) par(mfrow=c(3,2),bg="lightyellow") contour(myMvd1,dmvdc,xlim=c(0,9),ylim=c(0,11),col="blue", main="t-copula",xlab="duration",ylab="severity") points(D,S,col="red") contour(myMvd2,dmvdc,xlim=c(0,9),ylim=c(0,11),col="blue", main="gumbel-copula",xlab="duration",ylab="severity") points(D,S,col="red") #----------------Ud, Ui--------- myCop.t2<-tCopula(dim=2,param=0.8135,df=4.0545) myCop.g2<-gumbelCopula(2.2488) myMvd1<-mvdc(copula=myCop.t2,margins=c("weibull","gamma"), paramMargins=list(list(shape=1.9456,scale=3.616), list(shape=113.6474,rate=97.206))) myMvd2<-mvdc(copula=myCop.g2,margins=c("weibull","gamma"), paramMargins=list(list(shape=1.9456,scale=3.616), list(shape=113.6474,rate=97.206))) contour(myMvd1,dmvdc,xlim=c(0,9),ylim=c(0.75,1.6),col="blue", main="t-copula",xlab="duration",ylab="Intensity") points(D,I,col="red") contour(myMvd2,dmvdc,xlim=c(0,9),ylim=c(0.75,1.6),col="blue", main="gumbel-copula",xlab="duration",ylab="Intensity") points(D,I,col="red") #-------------------Ui, Us---------- myCop.t3<-tCopula(dim=2,param=0.8497,df=2.7614) myCop.g3<-gumbelCopula(2.5991) myMvd1<-mvdc(copula=myCop.t3,margins=c("gamma","gamma"), paramMargins=list(list(shape=113.6474,rate=97.206), list(shape=2.5196,rate=0.6482))) myMvd2<-mvdc(copula=myCop.g3,margins=c("gamma","gamma"), paramMargins=list(list(shape=113.6474,rate=97.206), list(shape=2.5196,rate=0.6482))) contour(myMvd1,dmvdc,xlim=c(0.75,1.6),ylim=c(0,11),col="blue", main="t-copula",xlab="Intensity",ylab="severity") points(I,S,col="red") contour(myMvd2,dmvdc,xlim=c(0.75,1.6),ylim=c(0,11),col="blue", main="gumbel-copula",xlab="Intensity",ylab="severity") points(I,S,col="red") ## calculate the loglikelihood for t-copula and gumbel copula data<-cbind(Ud,Us) loglikCopula(c(0.9955,1.4468e7),data,myCop.t1) # 108.0883 loglikCopula(c(0.9955,10),data,myCop.ts) loglikCopula(13.724,data,myCop.g1) #98.57883 gofCopula(myCop.t1,data,N=200,method="mpl") data<-cbind(Ud,Ui) loglikCopula(c(0.8135,4.0545),data,myCop.t2) # 24.79476 loglikCopula(2.2488,data,myCop.g2) #20.72553 data<-cbind(Ui,Us) loglikCopula(c(0.8497,2.7614),data,myCop.t3) #30.96441 loglikCopula(2.5991,data,myCop.g3) #26.80628 ############## Trivariate t copula ##### myCop.t<-tCopula(dim=3,param=c(0.9948,0.7561,0.8074),dispstr="un",df=7.2455) obs<-cbind(Ud,Us,Ui) pse<-apply(obs,2,rank)/47 #psudo-observation fit.ml<-fitCopula(myCop.t,pse,method="ml",start=c(0,0,0,10)) fit.ml t.copula <- tCopula(rep(0, 3), dim = 3, dispstr = "un", df.fixed=TRUE) gofCopula(t.copula,pse) ##-----conditional probability-----### d=1 s<-seq(0,10,by=0.5) F<-pweibull(d,shape=1.9456,scale=3.616) Fd<-rep(F,21) Fs<-pgamma(s,shape=2.5196,rate=0.6482) data=cbind(Fd,Fs) C<-pcopula(myCop.t1,cbind(Fd,Fs)) P<-(Fs-C)/(1-Fd) plot(s,P,pch=0,xlab="severity",ylab="probability") lines(s,P) #---- d=2 s<-seq(0,10,by=0.5) F<-pweibull(d,shape=1.9456,scale=3.616) Fd<-rep(F,21) Fs<-pgamma(s,shape=2.5196,rate=0.6482) data=cbind(Fd,Fs) C<-pcopula(myCop.t1,cbind(Fd,Fs)) P<-(Fs-C)/(1-Fd) points(s,P,pch=4) lines(s,P) #--- d=3 s<-seq(0,10,by=0.5) F<-pweibull(d,shape=1.9456,scale=3.616) Fd<-rep(F,21) Fs<-pgamma(s,shape=2.5196,rate=0.6482) data=cbind(Fd,Fs) C<-pcopula(myCop.t1,cbind(Fd,Fs)) P<-(Fs-C)/(1-Fd) points(s,P,pch=5) lines(s,P) #--- d=4 s<-seq(0,10,by=0.5) F<-pweibull(d,shape=1.9456,scale=3.616) Fd<-rep(F,21) Fs<-pgamma(s,shape=2.5196,rate=0.6482) data=cbind(Fd,Fs) C<-pcopula(myCop.t1,cbind(Fd,Fs)) P<-(Fs-C)/(1-Fd) points(s,P,pch=7) lines(s,P) #--- d=5 s<-seq(0,10,by=0.5) F<-pweibull(d,shape=1.9456,scale=3.616) Fd<-rep(F,21) Fs<-pgamma(s,shape=2.5196,rate=0.6482) data=cbind(Fd,Fs) C<-pcopula(myCop.t1,cbind(Fd,Fs)) P<-(Fs-C)/(1-Fd) points(s,P,pch=11) lines(s,P) #--- d=6 s<-seq(0,10,by=0.5) F<-pweibull(d,shape=1.9456,scale=3.616) Fd<-rep(F,21) Fs<-pgamma(s,shape=2.5196,rate=0.6482) data=cbind(Fd,Fs) C<-pcopula(myCop.t1,cbind(Fd,Fs)) P<-(Fs-C)/(1-Fd) points(s,P,pch=13) lines(s,P) #--- d=7 s<-seq(0,10,by=0.5) F<-pweibull(d,shape=1.9456,scale=3.616) Fd<-rep(F,21) Fs<-pgamma(s,shape=2.5196,rate=0.6482) data=cbind(Fd,Fs) C<-pcopula(myCop.t1,cbind(Fd,Fs)) P<-(Fs-C)/(1-Fd) points(s,P,pch=16) lines(s,P) #--- d=8 s<-seq(0,10,by=0.5) F<-pweibull(d,shape=1.9456,scale=3.616) Fd<-rep(F,21) Fs<-pgamma(s,shape=2.5196,rate=0.6482) data=cbind(Fd,Fs) C<-pcopula(myCop.t1,cbind(Fd,Fs)) P<-(Fs-C)/(1-Fd) points(s,P,pch=17) lines(s,P) legend(0.4,0.95,c("d'=1 month","d'=2 month","d'=3 month", "d'=4 month","d'=5 month","d'=6 month","d'=7 month","d'=8 month"), pch=c(0,4,5,7,11,13,16,17),bg='gray90') #-----measure goodness of fit------ ### goodness of fit t-copula--------- eu<-cbind((rank(D)-0.5)/46,(rank(S)-0.5)/46,(rank(I)-0.5)/46) C_empirical<-pcopula(myCop.t,eu) data<-cbind(Ui,Us,Ui) C_theoretical<-pcopula(myCop.t,data) par(mfrow=c(1,2),bg="lightcyan") qqplot(C_theoretical,C_empirical,main=" Goodness of fit t-copula") x<-seq(0,1,by=0.5) y<-seq(0,1,by=0.5) lines(x,y,col="red") ### goodness of fit gumbel-copula--- eu1<-cbind((rank(D)-0.5)/46,(rank(S)-0.5)/46) ec2<-pcopula(gumbelCopula(13.72),eu1) eu<-cbind(ec2,(rank(I)-0.5)/46) C_empirical<-pcopula(gumbelCopula(2.44),eu) c2<-pcopula(gumbelCopula(13.72),cbind(Ud,Us)) C_theoretical<-pcopula(gumbelCopula(2.44),cbind(c2,Ui)) qqplot(C_theoretical,C_empirical,main=" Goodness of fit Gumbel-copula") x<-seq(0,1,by=0.5) y<-seq(0,1,by=0.5) lines(x,y,col="red")