% Matlab code for Copula % data import dta=xlsread('drought events.xlsx'); D=dta(2:end,2); S=dta(2:end,3); I=dta(2:end,4); arrT=dta(2:end-1,6); %interarrival time % fit duration to weibull Ud=wblcdf(D,3.616,1.9456); % fit severity to gamma Us=gamcdf(S,2.5196,1.5428); % fit intensity to gamma Ui=gamcdf(I,113.7076,0.0103); %-----Parameter Estimation for bivariate copulas---------- % fit t-copula [rho1,nu1]=copulafit('t',[Ud Us]); [rho2,nu2]=copulafit('t',[Ud Ui]); [rho3,nu3]=copulafit('t',[Ui Us]); % fit gumbel param1=copulafit('Gumbel',[Ud Us]); param2=copulafit('Gumbel',[Ud Ui]); param3=copulafit('Gumbel',[Ui Us]); %---------Parameter Estimation for Trivariate copulas------ % fit t-copula [rho,nu]=copulafit('t',[Ud Us Ui]); % fit GH-copula theta1=copulafit('Gumbel',[Ud Us]); C2=copulacdf('Gumbel',[Ud Us],theta1); theta2=copulafit('Gumbel',[C2 Ui]); %------------------------ %%%------------Gumbel-----Bivariate Return period------ d=[0:0.1:8]; s=[0:0.1:12]; Fd=wblcdf(d,3.616,1.9456); Fs=gamcdf(s,2.5196,1.5428); EL=23.7333; for i=1:length(d) for j=1:length(s) Fds(i,j)=copulacdf('Gumbel',[Fd(i) Fs(j)],13); test(i,j)=1-Fd(i)-Fs(j)+Fds(i,j); Tds(i,j)=EL/test(i,j)/12;%year end end v=[3 5 10 20 50 100]; subplot(1,2,1) [C,h]=contour(s,d,Tds,v); clabel(C,h); xlabel('Severity'); ylabel('Duration (month)'); title('Return Periods T(D>=d and S>=s) from Gumbel') %------------ for i=1:length(d) for j=1:length(s) Fds(i,j)=copulacdf('Gumbel',[Fd(i) Fs(j)],13); test(i,j)=1-Fds(i,j); Tds(i,j)=EL/test(i,j)/12;%year end end subplot(1,2,2) [C,h]=contour(s,d,Tds,v); clabel(C,h); xlabel('Severity'); ylabel('Duration (month)'); title('Return Periods T(D>=d or S>=s) from Gumbel') %%%----------t---------Trivariate Return period---- % in order to compare single variable return period, using the same % quantile q1=0.604; q2=0.802; q3=0.901; q4=0.96; q5=0.98; q6=0.998; % compare return period=5 yr all R5=EL/(1-q1-q1-q1+copulacdf('t',[q1,q1],rho1,nu1)+... copulacdf('t',[q1,q1],rho2,nu2)+copulacdf('t',[q1,q1],rho3,nu3)... -copulacdf('t',[q1,q1,q1],rho,nu))/12; % 10 year R10=EL/(1-q2-q2-q2+copulacdf('t',[q2,q2],rho1,nu1)+... copulacdf('t',[q2,q2],rho2,nu2)+copulacdf('t',[q2,q2],rho3,nu3)... -copulacdf('t',[q2,q2,q2],rho,nu))/12; % 20 year R20=EL/(1-q3-q3-q3+copulacdf('t',[q3,q3],rho1,nu1)+... copulacdf('t',[q3,q3],rho2,nu2)+copulacdf('t',[q3,q3],rho3,nu3)... -copulacdf('t',[q3,q3,q3],rho,nu))/12; % 50 year R50=EL/(1-q4-q4-q4+copulacdf('t',[q4,q4],rho1,nu1)+... copulacdf('t',[q4,q4],rho2,nu2)+copulacdf('t',[q4,q4],rho3,nu3)... -copulacdf('t',[q4,q4,q4],rho,nu))/12; % 100 year R100=EL/(1-q5-q5-q5+copulacdf('t',[q5,q5],rho1,nu1)+... copulacdf('t',[q5,q5],rho2,nu2)+copulacdf('t',[q5,q5],rho3,nu3)... -copulacdf('t',[q5,q5,q5],rho,nu))/12; % 1000 year R1000=EL/(1-q6-q6-q6+copulacdf('t',[q6,q6],rho1,nu1)+... copulacdf('t',[q6,q6],rho2,nu2)+copulacdf('t',[q6,q6],rho3,nu3)... -copulacdf('t',[q6,q6,q6],rho,nu))/12; %%%%%--------or r5=EL/(1-copulacdf('t',[q1,q1,q1],rho,nu))/12; r10=EL/(1-copulacdf('t',[q2,q2,q2],rho,nu))/12; r20=EL/(1-copulacdf('t',[q3,q3,q3],rho,nu))/12; r50=EL/(1-copulacdf('t',[q4,q4,q4],rho,nu))/12; r100=EL/(1-copulacdf('t',[q5,q5,q5],rho,nu))/12; r1000=EL/(1-copulacdf('t',[q6,q6,q6],rho,nu))/12; %-----------------------------------