(*:Mathematica Version: 4.0  *)

(* :Copyright: Enrique Castillo *)

(*:Context: Statistics`Extremes`  *)

(*:Title: Extremes Functions  *)

(*:Authors: Enrique Castillo *)

(*:Summary:
This package provides basic Extreme value routines
*)


(*:Limitations:

*)

(*:Sources:   
E. Castillo, A. S. Hadi, N Balakrishnan and J. M. Sarabia.
Extremes ....
*)

(*:Keywords: Extremes *)

(*:Requirements: None.  *)

(*:Package Version: 1.0 *)

(*:History:   Created January 4, 2004 by E. Castillo
*)

BeginPackage["Statistics`Extremes`"];

Needs["Graphics`Graphics`"];
Needs["Statistics`ContinuousDistributions`"];

Options[ToUnitFrechetMarginal]={Estimates->False,Gumbel->False};

Options[ProbabilityPaperLine]={Grids->Automatic,PlotSize->400,\
ProbabilityTicks->Automatic};
Options[ProbabilityPaperLine1]={PlotSize->400};
Options[ProbabilityPaperLine2]={PlotSize->400};

MatrixTensorProduct::usage=
"MatrixTensorProduct[A,B] gives a new matrix that is the tensor product of \
matrices A and B. The matrices A and B can be of any dimension."

DataNames::usage=
"The list of names of the data sets in the book."

ShowDataNames::usage=
"Show the data names and their corresponding numbers."

MaximumDataIndices::usage=
"List of the indices of data sets in which the interest is on the right tail \
(maxima)."

MinimumDataIndices::usage=
"List of the indices of data sets in which the interest is on the left tail \
(minima)."

TableAlphaValues::usage=
"A matrix that gives the values of the error type I for several sample \
sizes."

DataSet::usage=
"DataSet[dataset] returns the indicated dataset as a list. dataset can be a \
number or a name."

Distr::usage=
"List of available probability papers."

ShowProbabilityPapers::usage=
"ShowProbabilityPapers returns the list of available probability papers."

ProbabilityPaperLine::usage=
"ProbabilityPaperLine[Sample,paper,lambda,NpointsLine] plots the given Sample \
on the indicated probability paper for the parameter lambda (only for Frechet \
and Weibull) and a regression line using NpointsLine on the tail of interest \
(if NpointsLine = 0, all points are used). The parameter paper can be a \
number or a name. The PlotSize is an optional parameter. The option Grids \
allows giving the desired grid lines and marks for X and the option PlotSize \
permits fixing the size of the plot."
 
ProbabilityPaperLine1::usage=
"ProbabilityPaperLine1[Sample,paper,lambda,NpointsLine,PlotSize:600] plots \
the given Sample on the indicated probability paper for the parameter lambda \
(only for Frechet and Weibull) and a regression line using NpointsLine on the \
tail of interest (if NpointsLine = 0, all points are used). The parameter \
paper can be a number or a name. The PlotSize is an optional parameter."
 
 ProbabilityPaperLine2::usage=
"ProbabilityPaperLine2[Sample,paper,lambda,NpointsLine,PlotSize:600] plots \
the given Sample on the indicated probability paper for the parameter lambda \
(only for Frechet and Weibull) and a regression line using NpointsLine on the \
tail of interest (if NpointsLine = 0, all points are used). The parameter \
paper can be a number or a name. The PlotSize is an optional parameter."

RoundNumber::usage=
"RoundNumber[x,factor] rounds the number x with decimals as indicated by \
factor."

FindAlpha::usage=
"FindAlpha[n] finds the value of alpha associated with a sample size of n in \
TableAlphaValues."

InitialGumbelMaximaEstimate::usage=
"InitialGumbelMaximaEstimate[Sample] returns initial max-Gumbel estimates for \
Sample."

InitialGumbelMinimaEstimate::usage=
"InitialGumbelMinimaEstimate[Sample] returns initial min-Gumbel estimates for \
Sample."

FGEVM::usage=
"FGEVM[x,lambda,delta,kappa] gives the cdf of the maximal GEV distribution."

FGEVMInv::usage=
"FGEVMInv[p,lambda,delta,kappa] returns the quantile of the maximal GEV \
distribution."

FGEVm::usage=
"FGEVm[x,lambda,delta,kappa] returns the cdf of the minimal GEV \
distribution."

FGEVmInv::usage=
"FGEVmInv[p,lambda,delta,kappa] returns the quantile of the minimal GEV \
distribution."

fGEVM::usage=
"fGEVM[x,lambda,delta,kappa] gives the pdf of the maximal GEV distribution."

fGEVm::usage=
"fGEVm[x,lambda,delta,kappa] gives the pdf of the minimal GEV distribution."

FGPDM::usage=
"FGPDM[x,alpha,kappa] returns the cdf of the maximal GPD distribution."

fGPDM::usage=
"fGPDM[x,alpha,kappa] returns the pdf of the maximal GPD distribution."

FGPDMInv::usage=
"FGPDMInv[p,alpha,kappa] returns the quantile of the maximal GPD \
distribution."

FGPDm::usage=
"FGPDm[x,alpha,kappa] gives the cdf of the minimal GPD distribution."

fGPDm::usage=
"fGPDm[x,alpha,kappa] gives the pdf of the minimal GPD distribution."

FGPDmInv::usage=
"FGPDmInv[p,alpha,kappa] returns the quantile of the minimal GPD \
distribution."

LKGEVM::usage=
"LKGEVM[Sample,lambda,delta,kappa,minVal,maxVal] returns the log likelihhod \
of the maximal GEV distribution for Sample."

GM::usage=
"GM[x,lambda,delta] gives the cdf of the maximal Gumbel distribution."

gM::usage=
"gM[x,lambda,delta] gives the pdf of the maximal Gumbel distribution."

Gm::usage=
"Gm[x,lambda,delta] gives the cdf of the minimal Gumbel distribution."

gm::usage=
"gm[x,lambda,delta] gives the pdf of the minimal Gumbel distribution."

LKGM::usage=
"LKGM[Sample,lambda,delta] gives the log likelihhod of the maximal Gumbel \
distribution for Sample."

LKGm::usage=
"LKGm[Sample,lambda,delta] gives the log likelihhod of the minimal Gumbel \
distribution for Sample."

HessGEVM::usage=
"HessGEVM[Sample,parameters] returns the covariance matrix of the ML \
estimates of the maximal GEV distribution."

EstimateGumbelML::usage=
"EstimateGumbelML[Sample,Tail,p] returns the Sample ML estimates of the \
Gumbel model for maxima (Tail = 1) or minima (Tail = 0). The p is an optional \
parameter for printing (Tail = 1) or not (Tail = 0)."

EstimateGEVML::usage=
"EstimateGEVML[Sample,Tail,p] returns the Sample ML estimates of the GEV for \
maxima (Tail = 1) or minima (Tail = 0). The p is an optional parameter for \
printing (Tail = 1) or not (Tail = 0)."

EstimateGEVMLThreshold::usage=
"EstimateGEVMLThreshold[TruncatedSample_,Tail_,u_,n_,p_:1] returns the Sample \
ML estimates of the GEV for maxima (Tail = 1) or minima (Tail = 0) of a \
truncated sample of size n at u. The p is an optional parameter for printing \
(Tail = 1) or not (Tail = 0)."

HessGPDM::usage=
"HessGPDM[Sample,parameters] returns the covariance matrix of the ML \
estimates of the maximal GPD distribution."

LKGPDM::usage=
"LKGPDM[Sample,alpha,kappa,maxVal] gives the log likelihhod of the maximal \
GPD distribution for Sample."

EstimateGPDML::usage=
"EstimateGPDML[Sample,Tail,p] returns the Sample ML estimates of the GPD for \
maxima (Tail = 1) or minima (Tail = 0). The p is an optional parameter for \
avoiding printing (p = 0)"

PPPlotGEVBands::usage=
"PPPlotGEVBands[Sample,lambda,delta,kappa,Tail,Title] plots the P-P plot \
associated with Sample and its ML GEV estimate corresponding to the indicated \
Tail: maxima (Tail = 1) or minima (Tail = 0)."

PPPlotGEV::usage=
"PPPlotGEV[Sample,lambda,delta,kappa,Tail,Title] plots the P-P plot \
associated with Sample and its ML GEV estimate corresponding to the indicated \
Tail: maxima (Tail = 1) or minima (Tail = 0)."

QQPlotGEV::usage=
"QQPlotGEV[Sample,lambda,delta,kappa,Tail,Title] plots the Q-Q plot \
associated with Sample and its ML GEV estimate corresponding to the indicated \
Tail: maxima (Tail = 1) or minima (Tail = 0)."

PPPlotGPDBands::usage=
"PPPlotGPDBands[Sample,alpha,kappa,Tail,Title] plots the P-P plot associated \
with Sample and its ML GPD estimate corresponding to the indicated Tail."

PPPlotGPD::usage=
"PPPlotGPD[Sample,alpha,kappa,Tail,Title] plots the P-P plot associated with \
Sample and its ML GPD estimate corresponding to the indicated Tail: maxima \
(Tail = 1) or minima (Tail = 0)."

QQPlotGPD::usage=
"QQPlotGPD[Sample,alpha,kappa,Tail,threshold,Title] plots the Q-Q plot \
associated with Sample and its ML GPD estimate corresponding to the indicated \
Tail: maxima (Tail = 1) or minima (Tail = 0). Threshold value is added to GPD \
sample."

PlotGEVM::usage=
"PlotGEVM[lambda,delta,kappa,p1,p2] plots the pdf of the maximal GEV \
distribution with the given parameters. The p1 and p2 are optional quantiles \
that define the range of the plot."

PlotGEVm::usage=
"PlotGEVm[lambda,delta,kappa,p1,p2] plots the pdf of the minimal GEV \
distribution with the given parameters. The p1 and p2 are optional quantiles \
that define the range of the plot."

PlotGPDM::usage=
"PlotGPDM[alpha,kappa,p1,p2] plots the pdf of the maximal GPD distribution \
with the given parameters. The p1 and p2 are optional quantiles that define \
the range of the plot."

PlotGPDm::usage=
"PlotGPDm[alpha,kappa,p1,p2] plots the pdf of the minimal GPD distribution \
with the given parameters. The p1 and p2 are optional quantiles that define \
the range of the plot."

CompleteGEVAnalysis::usage=
"CompleteGEVAnalysis[Sample,Tail,Title,PW] performs a complete extreme \
analysis of Sample with ML estimation probability plots and P-P and Q-Q \
Validation. The PW is an optional parameter that allows (PW = 1) calculating \
the Prescott-Walden asymtotic covariance matrix."

Curvature::usage=
"Curvature[Sample,Tail] calculates the given sample S statistic for the \
curvature method and the given Tail: maxima (Tail = 1) or minima (Tail = 0)."

CurvatureP::usage=
"CurvatureP[n_,NSimulations_,Tail_,cval_] calculates the cdf at the point \
cval for the curvature statistic and the corresponding p-value for a sample \
size n and with NSimulations."
 
fOrderStatistic::usage=
"fOrderStatistic[x,r,n] returns the pdf of the r-th Order statistic of a \
sample of size n from a parent with pdf f and cdf F (f and F are global \
functions defined by the user)."

FOrderStatistic::usage=
"FOrderStatistic[x,r,n] returns the cdf of the r-th Order statistic of a \
sample of size n from a parent with cdf F (F is a global function defined by \
the user)."

fTwoOrderStatistic::usage=
"FOrderStatistic[x,y,n,r,s] returns the pdf of the r-th and s-th Order \
statistics of a sample of size n from a parent with cdf F (F is a global \
function defined by the user)."

PlotOrderStatistic::usage=
"PlotOrderStatistic[f[x],F[x],r,n,xmin,xmax] plots the cdf of the r-th Order \
statistic of a sample of size n from a parent with pdf f[x] and cdf F[x] in \
the interval (xmin,xmax)."

PlotTwoOrderStatistics::usage=
"PlotTwoOrderStatistics[f[x],F[x],n,r,s,xmin,xmax,ymin,ymax] plots the cdf of \
the r-th and s-th Order statistics of a sample of size n from a parent with \
pdf f[x] and cdf F[x] in the region defined by (xmin,xmax) and (ymin,ymax)."

RangeEstimateGPDML::usage=
"RangeEstimateGPDML[Sample,u0,u1,m,Tail] plots the values of the GPD alpha \
and kappa parameter ML estimates versus u in the range (u0,u1) using m \
u-values equally spaced and the data in the indicated tail."

ASAEGPD::usage=
"ASAEGPD[Sample,lambda,kappa,Tail] calculates the average scaled relative \
error fro the GPD."

RangeEstimateGPDCH::usage=
"RangeEstimateGPDCH[Sample,umin,umax,m,Tail] plots the value of the GPD alpha \
parameter EPM estimate versus u in the range (umin,umax) using m equally \
spaced u-values."

ToUnitFrechetMarginal::usage=
"ToUnitFrechetMarginal[Sample,Tail,Gumbel] transforms the Sample to a Unit \
Frechet marginal sample. If the Gumbel option ios True, it uses the \
corresponding Gumbel model. If the Estimates option is selected returns the \
GEV parameters."

EstimateLogisticML::usage=
"EstimateLogisticML[Sample] estimates the alpha parameter of the logistic \
bivariate model by the maximum likelihood method."

EstimateLogisticLS::usage=
"EstimateLogisticLS[Sample] estimates the alpha parameter of the logistic \
bivariate model by the least squares method."

PointProcMLLogistic::usage=
"PointProcMLLogistic[Sample,p:1] estimates the alpha parameter of the \
logistic bivariate model by the point process ML method."

PointProcMLDirichlet::usage=
"PointProcMLDirichlet[Sample,p:1] estimates the alpha parameter of the \
Dirichlet bivariate model by the point process ML method."

PointProcJointMLLogisticGumbel::usage=
"PointProcJointMLLogisticGumbel[Sample,p] gives the joint estimate of all \
parameters of the Logistic bivariate model by the point process ML method. \
Optionally one can give the initial estimates (they can be obtained by the \
previous methods LS or ML). It is only for Gumbel marginals."

PointProcJointMLLogisticGEV::usage=
"PointProcJointMLLogisticGEV[Sample,p] gives the joint estimate of all \
parameters of the Logistic bivariate model by the point process ML method. \
Optionally one can give the initial estimates (they can be obtained by the \
previous methods LS or ML)."

PlotLogistic::usage=
"PlotLogistic[Sample,lambda1,delta1,lambda2,delta2,alpha,Size] plots the \
sample and the contours of the Logistic model with the given parameters. Size \
is an optional plot size."

 EstimateGEVDirichletML::usage=
"EstimateGEVDirichletML[Sample,p:1] estimates the complete GEVDirichlet \
model. The p is an optional parameter to print partial results."
 
PlotDirichlet::usage=
"PlotDirichlet[Sample,lambda11,delta11,lambda21,delta21,alpha1,beta1] plots \
the sample and the contours of the Dirichlet model with the given \
parameters."

PlotHistogramGEV::usage=
"PlotHistogramGEV[Sample, lambda, delta, kappa, d:8]  Plots a histogram of \
the Sample and the density of the GEV model with the given parameters."
     
EstimateGEVEPM::usage=
"EstimateGEVEPM[Sample,Tail,Option] gives the GEV estimates using the EPM. If \
Option = 1 it uses the trimmed mean; otherwise, the median estimate is used."

EstimateGEVEPM1::usage=
"EstimateGEVEPM1[Sample,Tail,Option] gives the GEV estimates using the EPM. \
If Option = 1 it uses the trimmed mean; otherwise, the median estimate is \
used."

EstimateGEVLS::usage=
"EstimateGPDLS[Sample,Tail] returns the Sample Least Squares estimates of the \
GEV for maxima (Tail = 1) or minima (Tail = 0)."

EstimateGEVEPMThreshold::usage=
"EstimateGEVEPMThreshold[TruncatedSample,Tail,n] gives the GEV estimates \
using the EPM with the truncated sample TruncatedSample. The n is the total \
sample size."      
      
EstimateGEVPWM::usage=
"EstimateGEVPWM[Sample,Tail] returns the PWM estimates of the GEV."

PlotPDFOrderStatistic1::usage=
"PlotPDFOrderStatistic1[f[x],F[x],x,r,n,xmin,xmax,p:] plots the pdf of the \
r-th order statistics or the list of order statistics in r of a sample of \
size n from a parent with pdf f[x] and cdf F[x] in the range (xmin,xmax) and \
p is an optional parameter to eliminate the plot."

PlotPDFOrderStatistic::usage=
"PlotPDFOrderStatistic[F[x],x,r,n,xmin,xmax,p:] plots the pdf of the r-th \
order statistics or the list of order statistics in r of a sample of size n \
from a parent with cdf F[x] in the range (xmin,xmax) and p is an optional \
parameter to avoid the plot (p = 0)."

PlotCDFOrderStatistic::usage=
"PlotCDFOrderStatistic[F[x],x,r,n,xmin,xmax,p] plots the cdf of the r-th \
order statistics or the list of order statistics in r of a sample of size n \
from a parent with cdf F[x] in the range (xmin,xmax). The p is an optional \
parameter to avoid the plot."
     
PlotPDFTwoOrderStatistics1::usage=
"PlotPDFTwoOrderStatistics1[f[x],F[x],var,r,s,n,xmin,xmax,ymin,ymax,p] plots \
the pdf of the r-th and s-th order statistics of a sample of size n from a \
parent with pdf f[x] and cdf F[x] in the range (xmin,xmax) and (ymin,ymax) \
and p is an optional parameter to avoid the plot."

PlotPDFTwoOrderStatistics::usage=
"PlotPDFTwoOrderStatistics[F[x],var,r,s,n,xmin,xmax,ymin,ymax,p] plots the \
pdf of the r-th and s-th order statistics of a sample of size n from a parent \
with cdf F[x] in the range (xmin,xmax) and (ymin,ymax) and p is an optional \
parameter to avoid the plot (p = 0)."
     
EstimateGPDCH1::usage=
"EstimateGPDCH1[Sample,Tail,Option,p] returns the Sample EPM estimates of the \
GPD for maxima (Tail = 1) or minima (Tail = 0). The p is an optional \
parameter for printing (Tail = 0) or not (Tail = 0). If Option is 1 the \
trimmed mean is used; otherwise, the median is used. A reduced number of \
elemental percentiles is used."

EstimateGPDCH::usage=
"EstimateGPDCH[Sample,Tail,Option,p] returns the Sample EPM estimates of the \
GPD for maxima (Tail=1) or minima (Tail=0). p is an optional parameter for \
printing (Tail = 0) or not (Tail = 0). If Option is 1 the trimmed mean is \
used; otherwise, the median is used. All elemental percentiles are used."

EstimateGPDPWM::usage=
"EstimateGPDPWM[Sample,Tail] returns the PWM estimates of the GPD."

EstimateGPDLS::usage=
"EstimateGPDLS[Sample,Tail] returns the Sample Least Squares estimates of the \
GPD for maxima (Tail = 1) or minima (Tail = 0)."

MakeBivariateSample::usage=
"MakeBivariateSample[MarginalX,MarginalY] builds the bivariate sample from \
the marginals."
     
PlotBivariateSample::usage=
"PlotBivariateSample[Sample,Size] plots a bivariate sample at an optional \
size."
     
DomainofAttraction::usage=
"DomainofAttraction[xp,var] determines the value of kappa for minima and \
maxima given the percentile function xp[var] and the variable var in that \
function."
     
PlotGEVLogistic::usage=
"PlotGEVLogistic[Sample,lambda1,delta1,kappa1,lambda2,delta2,kappa2,alpha,\
Size] plots the Logistic model including the six parameters to obtain Unit \
Frechet marginals. The ImageSize is an optional parameter for the Plot size."

PlotGEVDirichlet::usage=
"PlotGEVDirichlet[Sample,lambda1,delta1,kappa1,lambda2,delta2,kappa2,alpha,\
beta,Size] plots the Dirichlet model including the six parameters to obtain \
Unit Frechet marginals.  The ImageSize is an optional parameter for the Plot \
size."

ASAEGEV::usage=
"ASAEGEV[Sample,lambda,delta,kappa,Tail] calculates the average scaled \
absoluted error of a GEV estimate."
    
(*************************************************************************)

Begin["`Private`"]

DataNames={"Wind","Flood","Wave","Telephone1","Epicenter","Chain",\
"Insulation","Fatigue","Precipitation","Houmb","Macon","Hawkinsville","Women",\
"Men","Telephone2","Bilbao","Wind1","Wind2","Car1","Car2"};

ShowDataNames:=Do[Print[i,"   ",DataNames[[i]]],{i,1,Length[DataNames]}];

ShowProbabilityPapers:=Do[Print[i,"   ",Distr[[i]]],{i,1,Length[Distr]}];
      
MaximumDataIndices={1,2,3,10,11,12,13,14,16,17,18,19,20};

MinimumDataIndices={4,5,6,7,8,9,15};

TableAlphaValues={{10,0.0168},{20,0.0110},{40,0.0078},{50,0.0070},{100,0.0050}\
,
{200,0.0040},{500,0.0030},{1000,0.0026}};
    
DataSet[dataset_]:=Module[{CC,j},
If[StringQ[dataset],j=Position[DataNames,dataset][[1,1]],j=dataset];
CC[1]={22.64,22.80,23.75,24.01,24.04,24.24,24.74,25.45,25.55,25.66,25.99,
26.63,26.69,26.88,26.89,27.12,27.43,27.69,27.71,28.12,28.58,28.88,29.12,
29.45,29.48,30.18,31.31,31.55,31.57,32.54,32.98,33.83,33.86,34.64,35.21,
36.82,37.23,38.09,38.26,38.82,38.96,38.90,42.99,43.66,44.61,45.24,47.91,
54.75,69.40,98.16};
CC[2]={24.21,26.46,29.48,30.32,31.60,32.88,33.03,33.63,35.14,35.23,35.59,
35.89,35.95,36.07,36.49,36.50,37.13,37.48,38.01,38.21,38.53,38.91,39.26,
39.45,40.32,40.36,40.49,40.69,41.03,41.05,41.54,42.62,42.82,42.91,43.05,
43.31,43.34,43.42,43.65,43.87,44.71,45.04,45.58,46.00,48.29,48.76,49.28,
49.43,50.17,50.45,50.73,51.90,52.54,52.94,54.01,57.84,60.10,61.95,67.76,
75.70};
CC[3]={2.91,3.74,4.09,5.88,6.42,6.93,7.21,7.92,8.26,8.79,9.17,9.50,9.62,10.00,\

10.14,10.28,10.45,10.77,11.65,11.65,11.82,12.27,12.68,13.28,13.46,13.88,
13.98,14.32,14.38,14.46,14.86,15.03,15.30,16.07,16.23,17.36,18.68,18.72,
19.44,20.09,21.06,21.13,21.53,21.80,23.15,24.75,25.45,28.13,29.95,37.19};
CC[4]=10000*{0.000060,0.000074,0.000112,0.000200,0.000221,0.000236,0.000285,
0.000298,0.000337,0.000374,0.000389,0.000416,0.000487,0.000559,
0.000632,0.000645,0.000813,0.000960,0.001100,0.001130,0.001170,
0.001280,0.001300,0.001350,0.001420,0.001560,0.001590,0.001720,
0.002210,0.002380,0.002480,0.002640,0.003390,0.003480,0.005260};
CC[5]={58.20,58.20,59.50,61.80,65.80,67.80,68.50,70.90,73.70,77.00,80.80,
83.70,84.30,89.00,97.60,98.30,99.60,101.40,105.10,105.80,106.70,119.10,
119.50,119.90,121.90,125.70,128.40,146.10,153.90,154.60,155.80,157.40,
157.70,163.70,172.70,173.90,174.20,175.10,176.00,178.70,179.10,179.50,
180.70,182.10,182.70,186.70,187.50,191.00,192.60,193.00,199.40,211.60,
212.10,216.80,222.90,227.30,229.40,234.50,236.80,238.90};
CC[6]={51.1,57.1,65.1,69.2,71.2,73.1,76.1,81.9,84.1,86.6,88.1,92.6,94.9,96.6,
97.0,101.4,103.4,103.8,105.2,119.1};
CC[7]={744,822,847,885,920,948,968,985,1010,1018,1019,1028,1029,1031,1040,
1047,1071,1074,1097,1134,1147,1170,1174,1209,1251,1273,1320,1383,1388,
1462};
CC[8]={39611,44132,44209,45898,50139,54625,58970,64703,64950,66508,70208,
72098,75001,80393,81868,82202,82447,89268,90021,96136,96723,101610,
101833,106055,112833,119154,122366,134511,135220,136395,138378,153790,
184916,216370,240316}/1000000;
CC[9]={29.34,29.88,32.20,33.03,33.27,34.04,34.95,35.15,35.45,36.77,37.78,
38.35,38.37,39.14,39.52,40.00,40.47,40.48,41.05,41.15,41.75,42.62,43.05,
43.36,44.46,44.82,44.85,45.4,45.84,45.95,46.00,46.06,46.62,47.79,47.87,
48.13,49.06,49.42,49.63,52.13};
CC[10]={5.60,6.55,6.65,7.35,7.80,7.90,8.00,8.50,9.05,9.15,9.40,9.60,9.80,9.90,\

10.85,10.90,11.10,11.30,11.30,11.55,11.75,12.85,12.90,13.40};
CC[11]={28.8,8.5,44.8,51,4.8,19.1,47.8,25.4,14.3,31,66.2,37,48.6,28.3,21,72.5,\
28.3,
7.9,47.1,73.4,64.4,10.7,19.6,19,16.9,22.7,65.3,33.3,31,33.9,14.2,7.3,73.4,44.\
8,
50.2,40.4,57.6,32.6,24,84};
CC[12]={18.8,5.9,44.4,52,5.9,20.1,40.4,27,14.3,40,45.2,30,44,30.3,15.2,79.0,\
19.3,7.6,
42.4,70.5,50,12.2,16.2,19.9,17.4,13.5,61,25.8,33,37.9,13.3,6.9,57,41.6,46.8,
26.2,35.4,34.8,28.2,68};
CC[13]={101.50,101.69,102.12,102.32,102.54,102.72,102.78,102.92,103.14,103.31,\

103.38,103.41,103.46,103.53,103.56,103.56,103.77,103.83,103.86,103.94,
103.97,104.01,104.12,104.27,104.33,104.37,104.40,104.42,104.46,104.52,
104.71,104.85,104.87,105.01,105.01,105.02,105.03,105.19,105.32,105.45,
105.64,105.71,105.83,105.86,105.87,105.88,105.98,106.13,106.15,106.15,
106.52,107.49,107.89,107.90};
CC[14]={100.08,100.49,100.82,100.88,100.90,101.17,101.26,101.41,101.63,101.66,\

101.67,101.70,101.76,102.41,102.52,102.54,102.55,102.57,102.57,102.61,
102.63,102.69,102.78,102.88,102.94,103.00,103.06,103.15,103.17,103.24,
103.25,103.36,103.40,103.43,103.47,103.57,103.80,103.98,104.01,104.22,
104.65,104.88,104.92,105.00,105.12,105.12,105.12,105.48,105.55,105.72,
105.83,106.09,106.48,106.50};
CC[15]={1.34,0.14,0.33,1.68,1.86,1.31,0.83,0.33,2.20,0.62,3.20,1.38,0.96,0.28,\

0.44,0.59,0.25,0.51,1.61,1.85,0.47,0.41,1.46,0.09,2.18,0.07,0.02,0.64,
0.28,0.68,1.07,3.25,0.59,2.39,0.27,0.34,2.18,0.41,1.08,0.57,0.35,0.69,
0.25,0.57,1.90,0.56,0.09,0.28};
CC[16]={7.05,7.12,7.15,7.18,7.19,7.20,7.20,7.20,7.20,7.25,7.26,7.27,7.28,7.30,\

7.31,7.31,7.32,7.33,7.37,7.40,7.46,7.46,7.47,7.48,7.48,7.52,7.54,7.55,
7.55,7.58,7.59,7.59,7.61,7.63,7.65,7.66,7.66,7.67,7.67,7.68,7.69,7.72,
7.72,7.72,7.72,7.72,7.77,7.77,7.79,7.79,7.82,7.83,7.83,7.83,7.84,7.85,
7.85,7.88,7.88,7.90,7.90,7.91,7.93,7.93,7.93,7.94,7.95,7.95,7.97,7.97,
7.97,7.99,8.00,8.03,8.03,8.05,8.06,8.06,8.07,8.10,8.11,8.12,8.15,8.15,
8.15,8.18,8.18,8.18,8.19,8.20,8.21,8.23,8.23,8.30,8.30,8.31,8.31,8.32,
8.32,8.33,8.40,8.41,8.42,8.43,8.43,8.45,8.48,8.49,8.50,8.50,8.51,8.52,
8.53,8.54,8.56,8.58,8.59,8.59,8.60,8.65,8.69,8.71,8.72,8.74,8.74,8.74,
8.74,8.79,8.81,8.84,8.85,8.86,8.88,8.88,8.94,8.98,8.98,8.99,9.01,9.03,
9.06,9.12,9.16,9.17,9.17,9.18,9.18,9.18,9.21,9.22,9.23,9.24,9.27,9.29,
9.30,9.32,9.33,9.36,9.38,9.43,9.46,9.47,9.59,9.59,9.60,9.61,9.62,9.63,
9.66,9.74,9.75,9.78,9.79,9.79,9.80,9.84,9.85,9.89,9.90};
CC[17]={101.2,150.6,101.7,123.9,86.2,94.7,114.7,82.2,96.5,98.,70.1,100.9,108.\
2,101.2,
135.7,113.4,80.2,99.4,88.3,125.5,74.3,80.9,116.7,87.1,92.8,81.8,83.4,91.2,80.\
9,105.2,
90.8,83.4,102.,79.8,89.8,99.8,123.,96.2,91.7,111.2,108.6,129.6,76.3,
120.,104.9,110.3,84.6,98.9,85.3,99.8,70.8,81.8,77.1,123.6,102.1,141.,93.6,109.\
4,
86.2,90.2,98.2,84.8,119.1,99.6,113.5,87.7,111.1,106.5,78.9,97.3,97.5,79.9,
119.4,133.1,96.3,139.4,110.5,73.2,120.5,174.6,88.,86.7,100.1,101.2,109.3,73.0,\

107.6,71.6,133.7,107.5,134.8,87.9,114.7,82.5,91.6,81.7,136.8,85.,122.,83.4,82.\
0,
82.9,67.5,83.5,118.7,126.1,91.8,99.7,77.7,119.5,104.1,98.1,80.7,79.7,107.1,
84.8,88.5,132.9,88.,75.1,90.3,105.3,92.6,117.9,82.5,95.3,92.1,102.2,88.5,99.2,\

66.7,73.4,84.1,106.3,88.5,95.1,75.3,98.9,100.4,85.,107.9,101.7,86.7,81.9,89.8,\

130.4,117.9,88.2,81.5,92.9,92.8,104.5,116.3,117.9,92.5,101.8,100.8,80.,77.2,
80.1,97.4,98.5,95.,123.2,113.,130.5,83.,106.4,93.,58.1,89.1,91.,86.,98.3,125.\
3,
105.9,92.6,105.8,94.7,75.6,113.9,98.3,90.7,129.,86.1,173.,123.2,71.5,94.1,
74.6,105.1,87.9,70.5,74.6,89.2,128.2,109.,95.5,92.9,83.6};
CC[18]={96.8,154.9,115.,112.5,87.1,89.,114.4,78.6,90.7,92.,76.6,93.1,112.4,93.\
1,
152.1,120.5,86.3,89.3,91.1,136.6,96.,84.3,107.6,91.9,95.7,89.6,79.9,92.7,76.2,\
105.7,
93.8,82.3,112.2,79.2,108.8,92.1,112.1,89.1,71.8,105.8,105.4,113.4,76.7,115.5,\
102.2,
107.9,79.6,100.9,91.1,101.5,86.1,73.8,81.,105.5,96.7,142.9,96.8,111.7,83.4,90.\
7,98.,
85.9,112.8,92.4,101.9,87.2,121.,112.4,84.2,115.3,101.7,76.,130.1,127.8,102.0,
125.8,105.5,73.2,118.9,170.4,91.1,
89.3,93.5,110.1,98.1,68.9,118.2,78.3,162.4,109.3,106.5,96.8,108.8,86.9,
87.2,86.6,131.7,83.,119.4,86.,71.8,71.8,83.7,78.7,114.4,122.3,96.,94.5,
76.,113.7,109.1,98.9,78.4,84.,111.7,89.4,103.4,147.8,88.8,64.3,92.2,106.9,
80.7,114.5,84.9,97.1,83.2,96.,82.7,92.3,69.4,91.3,95.6,114.7,90.5,92.7,
80.4,92.9,106.1,89.2,110.2,95.9,83.7,86.4,80.7,130.7,110.2,84.2,83.3,122.,
84.7,81.2,111.2,122.5,97.2,86.1,98.4,91.3,73.4,84.3,102.3,100.9,91.8,
136.1,99.3,126.3,74.1,100.7,115.5,69.6,96.,77.5,86.2,103.4,127.1,87.3,
84.4,125.2,106.8,96.8,116.4,87.2,96.5,123.6,98.4,169.,120.,73.1,101.6,
65.9,101.6,79.6,81.4,73.2,90.6,112.1,105.3,94.2,101.4,73.9};
CC[19]={68.9,88.2,108.2,97.7,152.3,81.7,98.4,88.6,87.5,83.4,67.7,124.4,124.8,\
83.2,
114.1,109.3,165.8,85.6,122.9,93.9,118.9,84.9,75.1,106.,117.,127.7,88.2,80.9,
118.7,108.5,109.4,115.,98.,97.1,108.5,128.,102.6,88.4,129.2,81.2,102.9,75.,71.\
8,
105.6,111.8,114.3,141.3,164.5,132.7,123.6,102.3,115.5,76.2,80.5,121.3,65.1,
96.6,100.9,79.6,76.2,120.4,86.7,83.1,126.6,83.8,83.9,81.5,83.5,92.6,91.4,83.,
83.8,94.9,95.6,98.2,90.8,91.5,89.3,104.1,90.6,120.9,107.7,97.3,87.3,64.1,93.2,\

88.7,90.7,111.1,84.,81.1,81.5,89.4,114.6,138.3,134.3,119.,88.1,91.4,98.2,93.8,\

142.2,107.,81.8,99.1,84.4,82.,120.2,98.9,111.5,117.,95.1,112.,98.3,104.5,88.0,\

93.,108.,106.8,90.1,98.9,92.1,88.1,113.3,78.5,93.4,99.9,101.7,97.1,83.6,91.9,
109.2,138.2,92.4,98.5,85.8,81.3,118.8,121.3,133.8,107.1,78.9,81.,99.,80.4,105.\
7,
116.2,110.,71.7,85.3,84.8,79.6,95.4,72.4,115.1,90.4,97.1,100.7,85.6,78.2,90.1,\

108.6,90.4,81.1,102.3,101.3,83.6,99.5,95.2,127.5,90.4,114.2,93.8,103.5,79.5,
104.5,139.4,101.4,87.4,87.4,92.4,68.5,72.8,82.6,74.3,85.8,103.4,71.,91.9,102.\
2,
93.3,93.,127.3,106.7,75.3,101.2,93.8,73.9,130.8,131.3};
CC[20]={51.6,56.6,63.3,61.1,66.8,54.5,60.9,56.9,46.3,51.,59.5,64.8,65.7,56.1,\
63.4,
67.0,69.,60.2,64.5,57.4,62.6,59.3,53.3,61.2,65.3,64.8,61.6,59.,65.1,65.6,60.6,\
62.6,
61.3,62.1,61.1,67.8,54.9,54.4,67.4,63.2,60.,53.1,51.8,62.6,60.1,65.,68.9,68.2,\

66.1,64.7,57.1,63.3,57.9,55.7,66.9,42.,55.4,63.3,55.7,48.7,61.6,57.6,58.5,
64.,59.,49.5,52.4,57.3,59.4,56.7,53.4,55.,55.,64.4,56.2,65.,57.,56.6,59.5,62.\
7,
63.1,62.9,61.9,59.7,43.2,52.,48.1,58.2,62.3,64.2,58.8,61.8,61.7,64.3,66.8,
66.6,65.4,59.6,56.2,60.5,57.,66.7,55.9,57.8,55.,58.8,53.3,67.6,65.2,64.,62.7,
58.8,63.3,61.7,61.7,57.8,60.8,57.2,56.8,47.4,51.4,57.1,60.2,63.1,57.3,52.2,62.\
4,
63.4,60.8,54.4,57.5,65.,67.4,60.8,55.7,55.4,50.,59.9,65.2,67.7,61.6,64.2,55.9,\

63.5,54.3,62.7,65.7,61.8,55.9,57.3,49.7,56.2,54.4,48.1,64.7,61.,59.2,60.5,
57.3,53.1,55.8,62.1,64.4,56.,64.5,60.5,55.2,62.2,58.8,67.3,59.3,62.4,61.2,59.\
5,
55.5,63.4,67.4,56.1,61.9,64.8,54.7,53.2,50.3,56.9,59.3,52.1,60.1,42.6,68.0,
62.6,64.5,57.1,65.8,64.7,51.3,51.9,62.3,53.8,67.2,65.2};
Return[CC[j]]
];

Distr={"Normal","LogNormal","MaxWeibull","MaxGumbel","MaxFrechet",
    "MinWeibull","MinGumbel","MinFrechet"};
      
ProbabilityPaperLine1[Sample_,paper_,lambda_,NpointsLine_,options___]:=
Module[{XTran,YTran,PlotPointPosition,
Probabilities={0.00001,0.0001,0.001,0.005,0.01,0.02,0.05,0.1,0.2,0.3,
0.4,0.5,0.6,0.7,0.8,0.9,0.95,0.98,0.99,0.995,0.999,0.9999,0.99999},
      Probabilities0,Probabilities1,ReturnPeriods,i,i1,i2,iii,val, a,b,
      NpointsLine1=NpointsLine,SampleSize=Length[Sample],
      ProbPoints,Points,Sample1,Sample2,Y,meanX,meanY,DD1,DD2,DD3,DD8,DD9,SXX,\
SXY,Error=0,Tlog={2,3,5,6,8},
      XTransf={1,2,3,1,4,4,1,3},label={"X","X","X","X","X","X","X","X"},MX,bb,\

      Bandplus,Bandminus,aux5,aux6,u,u1,sigma,alpha=5/100,x1,x2,y1,y2,xmin,
      xmax,ymin,ymax,auxx,optionslist},
      optionslist={options};
      PlotSize=PlotSize/.optionslist/.Options[ProbabilityPaperLine1];
      If[StringQ[paper],paper1=Position[Distr,paper][[1,1]],paper1=paper];
    i1=-1;
    Do[If[SampleSize/(SampleSize+1)<Probabilities[[i]],i1=i1+1],{i,1,
        Length[Probabilities]}];
    Probabilities=Take[Probabilities,{i1+1,Length[Probabilities]-i1}];
    PTran[k_,x_]:=If[k<3,Quantile[NormalDistribution[0,1],x],
        If[k<6,-Log[-Log[x]],If[k<9,Log[-Log[1-x]],0]]];
    PTran1[k_,x_]:=If[k<3,CDF[NormalDistribution[0,1],x],
        If[k<6,Exp[-Exp[-x]],If[k<9,1-Exp[-Exp[x]],0]]];
    XTran[k_,x_]:=Which[k==1,x,
        k==2,x,k==3,1/(lambda-x),k==4,x,k==5,x-lambda,k==6,
        x-lambda,k==7,x,k==8,1/(lambda-x),True,0];
    PlotPointPosition[i_,SampleSize_]:=N[i/(SampleSize+1)];
    If[NpointsLine1==0,NpointsLine1=SampleSize];
    If[paper1==3 || paper1==8,
      label[[paper1]]=StringJoin["1/(",ToString[lambda],"-x)"]];
    If[lambda>0,
      If[paper1==5 || paper1==6,
        label[[paper1]]=StringJoin["x-",ToString[lambda]]]];
    If[lambda<0,
      If[paper1==5|| paper1==6,
        label[[paper1]]=StringJoin["x",ToString[lambda]]]];
    Print["Plotting a ",Distr[[paper1]]," probability paper"];
    Sample1=Map[XTran[paper1,#]&,Sample];
    Sample1=Sort[Sample1];
    If[paper1==3 || paper1==8,
      If[Sample1[[1]]<0,
        Print["****** Erroneous lambda. Minimum value of lambda must be: ",
              Max[Sample]," ******"];Error=1]];
    If[paper1==5 || paper1==6,
      If[Sample1[[1]]<0,
        Print["****** Erroneous lambda. Maximum value of lambda must be: ",
            Min[Sample]," ******"];Error=1;]];
    If[Abs[NpointsLine1]>SampleSize || Abs[NpointsLine1]<2,
         Print[
        "Error in the number of points to be used for the regression line it 
must be 2 <= n <= ",SampleSize];Error=1];
    If[Error!=1 ,
      If[NpointsLine1>0,i1=SampleSize-NpointsLine1+1;i2=SampleSize,i1=1;
        i2=-NpointsLine1];
      NpointsLine1=Abs[NpointsLine1];
      Probabilities0=Probabilities;
      Probabilities1=Map[{N[PTran[paper1,#]],ToString[#]}&,Probabilities0];
      ReturnPeriods=
        If[paper1<6,
          Table[{Probabilities1[[i,1]],
              Floor[100/(1-Probabilities0[[i]])+0.5]/100.0},{i,1,
              Length[Probabilities1]}],
          Table[{Probabilities1[[i,1]],
              Floor[100/Probabilities0[[i]]+0.5]/100.0},{i,1,
              Length[Probabilities1]}]];
      ProbPoints=
        N[Map[PTran[paper1,#]&,
            Table[PlotPointPosition[i,SampleSize],{i,1,SampleSize}]]];
      Points=Table[{Sample1[[i]],ProbPoints[[i]]},{i,1,SampleSize}];
      (* Regression line *)    
      If[Intersection[Tlog,{paper1}]!={},Sample2=Log[Sample1],
        Sample2=Sample1];
      (* Print["Sample1=",Sample1];*)
      (*Print["Sample2=",Sample2];*) 
         meanX=Sum[Sample2[[i]],{i,i1,i2}]/NpointsLine1;
      meanY=Sum[ProbPoints[[i]],{i,i1,i2}]/NpointsLine1;
      SXX=Sum[(Sample2[[i]]-meanX)^2,{i,i1,i2}]/NpointsLine1;
      SXY=
        Sum[(Sample2[[i]]-meanX)*(ProbPoints[[i]]-meanY),{i,i1,i2}]/
          NpointsLine1;
      b=SXY/SXX;a=meanY-b*meanX;
      If[a>0,Print["Regression line: v=",b,"u+",a],Print["Regression line: \
v=",b,"u",a]];
      (* Print[{{Sample2[[1]],a+b*Sample2[[1]]},{Sample2[[SampleSize]],
                a+b*Sample2[[SampleSize]]}}]; *)
      (* 
        Print["i1=",i1," i2=",i2," NpointsLine=",NpointsLine1]; *)  
      ExtremesLine={{Sample2[[1]],a+b*Sample2[[1]]},{Sample2[[SampleSize]],
            a+b*Sample2[[SampleSize]]}};
      ymax=Probabilities1[[Length[Probabilities1],1]];
      ymin=Probabilities1[[1,1]];
      If[ymax<ExtremesLine[[2,2]],auxx=(ymax-a)/b;
        ExtremesLine[[2]]={auxx,ymax}];
      If[Intersection[Tlog,{paper1}]!={},
        ExtremesLine[[2,1]]=Exp[ExtremesLine[[2,1]]]];
      If[ymin>ExtremesLine[[1,2]],
        auxx=(ymin-a)/b;ExtremesLine[[1]]={auxx,ymin}];
      If[Intersection[Tlog,{paper1}]!={},
        ExtremesLine[[1,1]]=Exp[ExtremesLine[[1,1]]]];
      MX=Table[{1,Sample2[[i]]},{i,i1,i2}];
      bb=Inverse[Transpose[MX].MX];
      sigma=
        Sqrt[Sum[(ProbPoints[[i]]-a-b*Sample2[[i]])^2,{i,i1,i2}]/(i2-i1)];
      (* Print["sigma=",sigma," bb=",bb]; *)
      Bandplus={};
      Bandminus={};
      xmin=Sample2[[1]];
      xmax=Sample2[[SampleSize]];
      ymin=Probabilities1[[1,1]];
      ymax=Probabilities1[[Length[Probabilities1],1]];
      Do[If[Intersection[Tlog,{paper1}]!={},u1=Exp[u],u1=u];
        iii=(SampleSize+1)*PTran1[paper1,a+b*u];
        (* Print["iii=",iii," 
PTran1[paper1,a+b*u]=",PTran1[paper1,a+b*u]];  *)
                If[iii<SampleSize+0.95 && iii>0.05,
          aux5=PTran[paper1,
              Quantile[BetaDistribution[iii,SampleSize-iii+1],1-alpha/2]]; 
          If[u>=xmin && u<=xmax && 
              aux5>=ymin && aux5 <=ymax,
            AppendTo[Bandplus,{u1,aux5}]];
          
          aux6=PTran[paper1,
              Quantile[BetaDistribution[iii,SampleSize-iii+1],alpha/2]];
          If[u>=xmin && u<=xmax && 
              aux6>=ymin &&aux6 <=ymax,
            AppendTo[Bandminus,{u1,aux6}]];]
        ,{u,Sample2[[i1]],Sample2[[i2]],(Sample2[[i2]]-Sample2[[i1]])/100}];
       (* Print["Bandplus=",Bandplus]; *)
       (* 
        Print["Bandminus=",Bandminus]; *)
      (* 
        Print["a=",a," b=",b];*)
      (* Print["val=",val];*)
      (* 
        Print["Points=",Points];*)
      (* 
        Print["Probabilities0=",Probabilities0];*)
      (* 
        Print["ProbPoints=",ProbPoints]; *)
      (* 
        Print["ReturnPeriods=",ReturnPeriods]; *)
      (* 
        Print["Probabilities1=",Probabilities1]; *)
      (* 
        Print["ExtremesLine=",ExtremesLine]; *)
      (* 
        Print["Range=",{{Points[[1,1]],Points[[1,2]]},{Points[[SampleSize,1]
],Points[[SampleSize,2]]}}];*)
      
      If[Intersection[Tlog,{paper1}]!={},
        DD8=LogLinearListPlot[Bandplus,       
            FrameTicks->{Automatic,Probabilities1,Automatic,
                ReturnPeriods},Frame->True,
            GridLines->{Automatic,PTran[paper1,Probabilities]},
            Prolog->AbsolutePointSize[4],      
            FrameLabel->{label[[paper1]],"Probability",Distr[[paper1]],
                "Return Period"},PlotJoined->True,
            PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
            DisplayFunction:>Identity];
        DD9=LogLinearListPlot[Bandminus,
FrameTicks->{Automatic,Probabilities1,Automatic,
                ReturnPeriods},Frame->True,
            GridLines->{Automatic,PTran[paper1,Probabilities]},
            Prolog->AbsolutePointSize[4], 
            FrameLabel->{label[[paper1]],"Probability",Distr[[paper1]],
                "Return Period"},PlotJoined->True,
            PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
            DisplayFunction:>Identity];
        DD1=LogLinearListPlot[ExtremesLine,     
            FrameTicks->{Automatic,Probabilities1,Automatic,
                ReturnPeriods},Frame->True,
            GridLines->{Automatic,PTran[paper1,Probabilities]},
            Prolog->AbsolutePointSize[4],           
            FrameLabel->{label[[paper1]],"Probability",Distr[[paper1]],
                "Return Period"},PlotJoined->True,
            PlotStyle->AbsoluteThickness[1],
            DisplayFunction:>Identity];
        DD2=
          LogLinearListPlot[Points,
            FrameTicks->{Automatic,Probabilities1,Automatic,
                ReturnPeriods},
            AxesOrigin->{Points[[1,1]],Points[[1,2]]},
            GridLines->{Automatic,PTran[paper1,Probabilities]},
            Prolog->AbsolutePointSize[4],            
            FrameLabel->{label[[paper1]],"Probability",Distr[[paper1]],
                "Return Period"},DisplayFunction:>Identity],
        DD8=
          ListPlot[Bandplus,
            FrameTicks->{Automatic,Probabilities1,Automatic,
                ReturnPeriods},
            AxesOrigin->{Points[[1,1]],Points[[1,2]]},
            Frame->True,
            GridLines->{Automatic,PTran[paper1,Probabilities]},
            Prolog->AbsolutePointSize[4],         
            FrameLabel->{"X","Probability",Distr[[paper1]],
                "Return Period"},PlotJoined->True,
            PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
            DisplayFunction:>Identity];
        DD9=
          ListPlot[Bandminus,
            FrameTicks->{Automatic,Probabilities1,Automatic,
                ReturnPeriods},
            AxesOrigin->{Points[[1,1]],Points[[1,2]]},
            Frame->True,
            GridLines->{Automatic,PTran[paper1,Probabilities]},
            Prolog->AbsolutePointSize[4],
            
            FrameLabel->{"X","Probability",Distr[[paper1]],
                "Return Period"},PlotJoined->True,
            PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
            DisplayFunction:>Identity];
        DD1=
          ListPlot[ExtremesLine,
            FrameTicks->{Automatic,Probabilities1,Automatic,
                ReturnPeriods},
            AxesOrigin->{Points[[1,1]],Points[[1,2]]},
            Frame->True,
            GridLines->{Automatic,PTran[paper1,Probabilities]},
            Prolog->AbsolutePointSize[4],
            
            FrameLabel->{"X","Probability",Distr[[paper1]],
                "Return Period"},PlotJoined->True,
            PlotStyle->AbsoluteThickness[1],
            DisplayFunction:>Identity];
        DD2=
          ListPlot[Points,
            PlotRange->{{Points[[1,1]],
                  Points[[1,2]]},{Points[[SampleSize,1]],
                  Points[[SampleSize,2]]}},
            FrameTicks->{Automatic,Probabilities1,Automatic,
                ReturnPeriods},
            AxesOrigin->{Points[[1,1]],Points[[1,2]]},
            GridLines->{Automatic,PTran[paper1,Probabilities]},
            Prolog->AbsolutePointSize[4],          
            FrameLabel->{"X","Probability",Distr[[paper1]],
                "Return Period"},DisplayFunction:>Identity]
        ],Return[]];
    DD3=Show[DD1,DD2,DD8,DD9,Frame->True,
        DisplayFunction:>$DisplayFunction,PlotRange->All,
        ImageSize->{PlotSize,Automatic}];
    Return[DD3]
    ];

RoundNumber[x_,factor_]:=Floor[factor*x+0.5]/factor;

FindAlpha[n_]:=Module[{i=1,TableAlphaValues},
    TableAlphaValues={{0,0.01},{10,0.0168},{20,0.0110},{40,0.0078},{50,
          0.0070},{100,0.0050},{200,0.0040},{500,0.0030},{1000,0.0026},{10^9,
          0.0026}};
    imax=Length[TableAlphaValues];
    While[TableAlphaValues[[i,1]]<n ,
      i=i+1];i=i-1;
    Print[i];
    If[n==TableAlphaValues[[i+1,1]],Return[TableAlphaValues[[i+1,2]]],
      Return[
        TableAlphaValues[[i,
              2]]+(TableAlphaValues[[i+1,2]]-
                TableAlphaValues[[i,2]])/(TableAlphaValues[[i+1,1]]-n)]]
    ];
    
InitialGumbelMaximaEstimate[Sample_]:=
  Module[{n=Length[Sample],i,Var,mean,delta,lambda},
    mean=Sum[Sample[[i]],{i,1,n}]/n;
    Var=Sum[(Sample[[i]]-mean)^2,{i,1,n}]/n;
    delta=Sqrt[Var*6/Pi^2];
    lambda=mean-0.57772*delta;
    Return[{lambda,delta}]
    ];
    
InitialGumbelMinimaEstimate[Sample_]:=
  Module[{n=Length[Sample],i,Var,mean,delta,lambda},
    mean=Sum[Sample[[i]],{i,1,n}]/n;
    Var=Sum[(Sample[[i]]-mean)^2,{i,1,n}]/n;
    delta=Sqrt[Var*6/Pi^2];
    lambda=mean+0.57772*delta;
    Return[{lambda,delta}]
    ];
   
FGEVM[x_,lambda_,delta_,kappa_]:=
    If[kappa>0.05,
      If[(x<=lambda+delta/kappa) ,Exp[-(1-kappa*(x-lambda)/delta)^(1/kappa)],
        1],If[kappa<-0.05,
        If[(x>=lambda+delta/kappa) ,
          Exp[-(1-kappa*(x-lambda)/delta)^(1/kappa)],0],
        Exp[-Exp[(lambda-x)/delta]]]];
FGEVMInv[p_,lambda_,delta_,kappa_]:=
    If[Abs[kappa]>0.01,lambda+delta*(1-(-Log[p])^kappa)/kappa,
      lambda-delta*Log[-Log[p]]];
FGEVm[x_,lambda_,delta_,kappa_]:=
    If[kappa>0.05,
      If[(x>=lambda-delta/kappa) ,
        1-Exp[-(1+kappa*(x-lambda)/delta)^(1/kappa)],0],
      If[kappa<-0.05,
        If[(x<=lambda-delta/kappa) ,
          1-Exp[-(1+kappa*(x-lambda)/delta)^(1/kappa)],1],
        1-Exp[-Exp[(x-lambda)/delta]]]];
FGEVmInv[p_,lambda_,delta_,kappa_]:=
    If[Abs[kappa]>0.01,lambda-delta*(1-(-Log[1-p])^kappa)/kappa,
      lambda+delta*Log[-Log[1-p]]];
      
fGEVM[x_,lambda_,delta_,kappa_]=
    If[Abs[kappa]>0.01,
      Exp[-(1-kappa*(x-lambda)/delta)^(1/kappa)] *(1-
                kappa*(x-lambda)/delta)^(1/kappa-1)/delta,
      Exp[-Exp[(lambda-x)/delta]+(lambda-x)/delta]/delta];

fGEVm[x_,lambda_,delta_,kappa_]:=
    If[Abs[kappa]>0.01,
      Exp[-(1+kappa*(x-lambda)/delta)^(1/kappa)] *(1+
                kappa*(x-lambda)/delta)^(1/kappa-1)/delta,
      Exp[-Exp[(x-lambda)/delta]+(x-lambda)/delta]/delta];
      
FGPDM[x_,alpha_,kappa_]:=
    If[x<0,0,If[Abs[kappa]>0.01,
        If[((kappa>0 &&  x<alpha/kappa) || (kappa<0)),
          1-(1-kappa*x/alpha)^(1/kappa),1],1-Exp[-x/alpha]]];
          
fGPDM[x_,alpha_,kappa_]:=
    If[x<0,0,If[Abs[kappa]>0.01,
        If[((kappa>0 &&  x<alpha/kappa) || (kappa<0)),(1-
                  kappa*x/alpha)^(1/kappa-1)/alpha,1],Exp[-x/alpha]/alpha]];
                  
FGPDMInv[p_,alpha_,kappa_]:=
    If[Abs[kappa]>0.01,alpha*(1-(1-p)^kappa)/kappa,-alpha*Log[1-p]];
    
FGPDm[x_,alpha_,kappa_]:=
    If[x>0,1,If[Abs[kappa]>0.01,
        If[((kappa>0 &&  x>-alpha/kappa) || (kappa<0)) ,(1+kappa*x/alpha)^(1/
                kappa),0],Exp[x/alpha]]];
                
fGPDm[x_,alpha_,kappa_]:=
    If[x>0,1,If[Abs[kappa]>0.01,
        If[((kappa>0 &&  x>-alpha/kappa) || (kappa<0)) ,(1+
                  kappa*x/alpha)^(1/kappa-1)/alpha,0],Exp[x/alpha]/alpha]];
                  
FGPDmInv[p_,alpha_,kappa_]:=
    If[Abs[kappa]>0.01,alpha*(p^kappa-1)/kappa,alpha*Log[p]];
      
      LKGEVM[Sample_,lambda_,delta_,kappa_,minVal_,maxVal_]:=
    Module[{z,n,aux1=-10^6,i},n=Length[Sample];
      If[((lambda+delta/kappa>=maxVal&&kappa>0&&
                delta>0)||(lambda+delta/kappa<=minVal&&kappa<0&&delta>0)),
        z=Table[-Log[Max[1.-kappa*(Sample[[i]]-lambda)/delta,0.000001]]/
              kappa,{i,1,n}];
        aux1=-n*Log[delta]-(1-kappa)*Sum[z[[i]],{i,1,n}]-
            Sum[Exp[-z[[i]]],{i,1,n}]];
      Return[aux1]];
      
GM[x_,lambda_,delta_]:=Exp[-Exp[(lambda-x)/delta]];

gM[x_,lambda_,delta_]=D[GM[x,lambda,delta],x];

Gm[x_,lambda_,delta_]:=1-Exp[-Exp[(x-lambda)/delta]];

gm[x_,lambda_,delta_]=D[Gm[x,lambda,delta],x];

LKGM[Sample_,lambda_,delta_]:=If[delta>0,-Length[Sample]*Log[delta]-
      Sum[Exp[(lambda-Sample[[i]])/delta]+(Sample[[i]]-lambda)/delta,{i,1,
          Length[Sample]}],-10^6];
          
LKGm[Sample_,lambda_,delta_]:=
    Sum[Log[gm[Sample[[i]],lambda,delta]],{i,1,Length[Sample]}];
    
HessGEVM[Sample_,parameters_]:=
  Module[{eps,ep=0.00001,i,j,n,x1,x2,x3,x4,m,minVal,maxVal},
    $MinPrecision=15;
	$MaxPrecision=15;
    maxVal=Max[Sample];
    minVal=Min[Sample];
    n=Length[parameters];
    eps=-ep*parameters;
    m=Table[0,{i,1,n},{j,1,n}];
    Do[
      x1=parameters;
      x1[[i]]=x1[[i]]+eps[[i]];
      x1[[j]]=x1[[j]]+eps[[j]];
      x2=parameters;
      x2[[i]]=x2[[i]]+eps[[i]];
      x2[[j]]=x2[[j]]-eps[[j]];
      x3=parameters;
      x3[[i]]=x3[[i]]-eps[[i]];
      x3[[j]]=x3[[j]]+eps[[j]];
      x4=parameters;
      x4[[i]]=x4[[i]]-eps[[i]];
      x4[[j]]=x4[[j]]-eps[[j]];
      m[[i,
            j]]=-(LKGEVM[Sample,x1[[1]],x1[[2]],x1[[3]],minVal,maxVal]-
                LKGEVM[Sample,x2[[1]],x2[[2]],x2[[3]],minVal,maxVal]-
                LKGEVM[Sample,x3[[1]],x3[[2]],x3[[3]],minVal,maxVal]+
                LKGEVM[Sample,x4[[1]],x4[[2]],x4[[3]],minVal,maxVal])/(4*
              eps[[i]]*eps[[j]]);
      ,{j,1,n},{i,1,n}];
    Return[Inverse[m]]
    ];

EstimateGumbelML[Sample_,Tail_,p_:1]:=Module[{lambdaG,deltaG,Sample1,sol},
    If[Tail==0,Sample1=-Sample,Sample1=Sample];
    {lambdaG,deltaG}=InitialGumbelMaximaEstimate[Sample1];
    sol=FindMinimum[-LKGM[Sample1,lambda,delta],{lambda,lambdaG,
          lambdaG*1.1},{delta,deltaG,deltaG*1.1},MaxIterations->100];
    {lambdaG,deltaG}={lambda,delta}/.sol[[2]];
    If[Tail==0,lambdaG=-lambdaG];
    If[p==1,
      Print["Likelihood=",-sol[[1]]," lambda=",lambdaG," delta=",deltaG]];
    Return[{sol[[1]],lambdaG,deltaG}]
    ];
    
EstimateGEVML[Sample_,Tail_,p_:1]:=
    Module[{maxVal,minVal,lambda0,delta0,kappa0,lambdaG,deltaG,kappaG,sol,
        sol1,Var,se,ii,likelihoodG,likelihood,Sample1,out},
      If[Tail==0,Sample1=-Sample,Sample1=Sample];
      maxVal=Max[Sample1];
      minVal=Min[Sample1];
      {lambda0,delta0}=InitialGumbelMaximaEstimate[Sample1];
      sol1=
            FindMinimum[-LKGM[Sample1,lambda,delta],{lambda,lambda0,
                lambda0*1.1},{delta,delta0,delta0*1.1},
              MaxIterations->100];
        likelihoodG=sol1[[1]];
        {lambdaG,deltaG}={lambda,delta}/.sol1[[2]];
      {lambda0,delta0,kappa0}=EstimateGEVEPM1[Sample1,1,0];
      sol=
        FindMinimum[-LKGEVM[Sample1,lambda,delta,kappa,minVal,maxVal],{lambda,\

            lambda0,lambda0*1.1},{delta,delta0,delta0*1.1},{kappa,kappa0,
            kappa0*1.1},MaxIterations->100];
      likelihood=sol[[1]];
      {lambda0,delta0,kappa0}={lambda,delta,kappa}/.sol[[2]];
      Var=HessGEVM[Sample1,{lambda0,delta0,kappa0}];
      se={};
      Do[AppendTo[se,Sqrt[Var[[ii,ii]]]],{ii,1,Length[Var]}];
      If[Tail==0,lambdaG=-lambdaG;lambda0=-lambda0;
        Var[[1,2]]=-Var[[1,2]];Var[[2,1]]=-Var[[2,1]];Var[[1,3]]=-Var[[1,3]];
        Var[[3,1]]=-Var[[3,1]];];
      If[p==1,
        Print["Gumbel likelihood and parameter estimates: Likelihood=",
          likelihoodG," lambda=",lambdaG," delta=",deltaG];
        Print["GEV likelihood and parameter estimates: Likelihood=",
          likelihood," lambda=",lambda0," delta=",delta0," kappa=",kappa0];
        Print["Var=",MatrixForm[Var]];
        Print["Standard Error=",se]];
      Return[{likelihoodG,{lambdaG,deltaG},likelihood,{lambda0,delta0,kappa0},\

          Var,se}]];

LKGEVMThreshold[Sample_,lambda_,delta_,kappa_,u_,minVal_,maxVal_]:=
    Module[{z,n,aux1,i},n=Length[Sample];
      If[((lambda+delta/kappa>maxVal&&kappa>0&&
                delta>0)||(lambda+delta/kappa<minVal&&kappa<0&&delta>0)),
        z=Table[-Log[1.-kappa*(Sample[[i]]-lambda)/delta]/kappa,{i,1,n}];
        aux1=-n*Log[delta]-(1-kappa)*Sum[z[[i]],{i,1,n}]-
            Sum[Exp[-z[[i]]],{i,1,n}]-
            n*Log[Max[0.000001,1-Exp[-(1-kappa*(u-lambda)/delta)^(1/kappa)]]],\

        aux1=-10^6];
      Return[aux1]];

EstimateGEVMLThreshold[TruncatedSample_,Tail_,u_,n_,p_:1]:=
    Module[{maxVal,minVal,lambda0,delta0,kappa0,lambdaG,deltaG,kappaG,sol,
        sol1,Var,se,ii,likelihoodG,likelihood,Sample1,out},
      If[Tail==0,Sample1=-TruncatedSample,Sample1=TruncatedSample];
      maxVal=Max[Sample1];
      minVal=Min[Sample1];
      {lambdaG,deltaG,kappaG}=EstimateGEVEPMThreshold[Sample1,1,n,0];
       If[p==1,Print["Initial Estimates : ",{lambdaG,deltaG,kappaG}]];
      sol=
        FindMinimum[-LKGEVMThreshold[Sample1,lambda,delta,kappa,u,minVal,
              maxVal],{lambda,lambdaG,lambdaG*1.1},{delta,deltaG,
            deltaG*1.1},{kappa,kappaG,kappaG*1.1},MaxIterations->100];
      likelihood=sol[[1]];
      {lambda0,delta0,kappa0}={lambda,delta,kappa}/.sol[[2]];
      If[Tail==0,lambdaG=-lambdaG;lambda0=-lambda0];
      If[p==1,
        Print["GEV likelihood and parameter estimates: Likelihood=",
          likelihood," lambda=",lambda0," delta=",delta0," kappa=",kappa0]];
      Return[{likelihood,{lambda0,delta0,kappa0}}]];

LKGPDM[Sample_,alpha_,kappa_,maxVal_]:=
    Module[{n,aux1=-10^6,i},n=Length[Sample];
      If[Abs[kappa]<0.01,aux1=-n*Log[alpha]-Sum[Sample[[i]],{i,1,n}]/alpha,
        If[((alpha/kappa>maxVal&&kappa>0)||kappa<0),
          aux1=-n*Log[alpha]+(1/kappa-1)*
                Sum[Log[1-kappa*Sample[[i]]/alpha],{i,1,n}]]];
      Return[aux1]];

EstimateGPDML::invalid="`1` contains negative values.";

EstimateGPDML[Sample_,Tail_,p_:1]:=Module[{maxVal,minVal,alpha0,kappa0,sol,\
Var,se,ii,
n=Length[Sample],Sample1},
If[Tail==0,Sample1=-Sample,Sample1=Sample];
If[Min[Sample1]<0,If[Tail==1,Print["Invalid sample because it contains 
negative numbers"],
Print["Invalid sample because it contains positive numbers"]];
Return[]];
maxVal=Max[Sample1];
minVal=Min[Sample1];
alpha0=Sum[Sample1[[i]],{i,1,n}]/n;
sol=FindMinimum[-LKGPDM[Sample1,alpha,kappa,maxVal],{alpha,alpha0,alpha0*1.1},\

{kappa,0.1,-0.2},MaxIterations->100];
{alpha0,kappa0}={alpha,kappa}/.sol[[2]];
Var=HessGPDM[Sample1,{alpha0,kappa0}];
se={};
Do[AppendTo[se,Sqrt[Var[[ii,ii]]]],{ii,1,Length[Var]}];
If[p==1,
Print["GPD likelihood and parameter estimates: Likelihood=",sol[[1]]," 
alpha=",alpha0," kappa=",kappa0];
Print["Var=",MatrixForm[Var]];
Print["Standard deviations =",se];];
Return[{sol[[1]],{alpha0,kappa0},Var,se}]];
   
 HESSGPDM1[Sample_,parameters_]:=
    Module[{n,aux1={{0,0},{0,0}},i,alpha,kappa,alpha1,kappa1},{alpha1,kappa1}=\

        parameters;
      n=Length[Sample];
      maxVal=Max[Sample];
      If[Abs[kappa1]<0.01,
        aux1={{D[-n*Log[alpha1]-Sum[Sample[[i]],{i,1,n}]/alpha1,alpha,
                  alpha]/.{alpha->alpha1,kappa->kappa1},
              D[-n*Log[alpha1]-Sum[Sample[[i]],{i,1,n}]/alpha1,alpha,
                  kappa]/.{alpha->alpha1,
                  kappa->kappa1}},{D[-n*Log[alpha1]-
                    Sum[Sample[[i]],{i,1,n}]/alpha1,kappa,
                  alpha]/.{alpha->alpha1,kappa->kappa1},
              D[-n*Log[alpha1]-Sum[Sample[[i]],{i,1,n}]/alpha1,kappa,
                  kappa]/.{alpha->alpha1,kappa->kappa1}}},
        If[(kappa1>0&&((alpha1/kappa1>maxVal&&kappa1>0)||kappa1<0)),
          aux1={{D[-n*Log[alpha]+(1/kappa-1)*
                        Sum[Log[1-kappa*Sample[[i]]/alpha],{i,1,n}],alpha,
                    alpha]/.{alpha->alpha1,kappa->kappa1},
                D[-n*Log[alpha]+(1/kappa-1)*
                        Sum[Log[1-kappa*Sample[[i]]/alpha],{i,1,n}],alpha,
                    kappa]/.{alpha->alpha1,
                    kappa->kappa1}},{D[-n*Log[alpha]+(1/kappa-1)*
                        Sum[Log[1-kappa*Sample[[i]]/alpha],{i,1,n}],kappa,
                    alpha]/.{alpha->alpha1,kappa->kappa1},
                D[-n*Log[alpha]+(1/kappa-1)*
                        Sum[Log[1-kappa*Sample[[i]]/alpha],{i,1,n}],kappa,
                    kappa]/.{alpha->alpha1,kappa->kappa1}}}]];
      Print["Hessian=",-aux1];
      Return[Inverse[-aux1]]];
      
HessGPDM[Sample_,parameters_]:=
  Module[{eps,ep=0.0001,i,j,n,x1,x2,x3,x4,m,maxVal},
maxVal=Max[Sample];  
    n=Length[parameters];
    eps=-ep*parameters;
    m=Table[0,{i,1,n},{j,1,n}];
    Do[x1=parameters;
      x1[[i]]=x1[[i]]+eps[[i]];
      x1[[j]]=x1[[j]]+eps[[j]];
      x2=parameters;
      x2[[i]]=x2[[i]]+eps[[i]];
      x2[[j]]=x2[[j]]-eps[[j]];
      x3=parameters;
      x3[[i]]=x3[[i]]-eps[[i]];
      x3[[j]]=x3[[j]]+eps[[j]];
      x4=parameters;
      x4[[i]]=x4[[i]]-eps[[i]];
      x4[[j]]=x4[[j]]-eps[[j]];
      m[[i,
            j]]=-(LKGPDM[Sample,x1[[1]],x1[[2]],maxVal]-
                LKGPDM[Sample,x2[[1]],x2[[2]],maxVal]-
                LKGPDM[Sample,x3[[1]],x3[[2]],maxVal]+
                LKGPDM[Sample,x4[[1]],x4[[2]],maxVal])/(4*eps[[i]]*
              eps[[j]]);,{j,1,n},{i,1,n}];
    Return[Inverse[m]]];
  
PPPlotGEVBands[Sample_,lambda_,delta_,kappa_,Tail_,Title_]:=
  Module[{aux,SS=Sort[Sample],n=Length[Sample],aux1,i,alpha1},
    alpha1=FindAlpha[Length[Sample]];
    aux=Table[{i/(n+1),
          If[Tail==1,FGEVM[SS[[i]],lambda,delta,kappa],
            FGEVm[SS[[i]],lambda,delta,kappa]]},{i,1,n}];
    aux1={{0,0}};
    Do[AppendTo[
        aux1,{i/(n+1),Quantile[BetaDistribution[i,n-i+1.],alpha1/2]}],{i,1,
        n}];
    AppendTo[aux1,{1,1}];
    D1=ListPlot[aux1,
        PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
        PlotJoined->True,DisplayFunction->Identity,PlotRange->All];
    aux1={{0,0}};
    Do[AppendTo[
        aux1,{i/(n+1),Quantile[BetaDistribution[i,n-i+1.],1-alpha1/2]}],{i,1,
        n}];
    AppendTo[aux1,{1,1}];
    D2=ListPlot[aux1,
        PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
        PlotJoined->True,DisplayFunction->Identity,PlotRange->All];
    DD1=ListPlot[aux,Prolog->AbsolutePointSize[5],Axes->False,
        DisplayFunction:>Identity];
    DD2=Graphics[Line[{{0,0},{1,1}}]];
    aux1=Show[D1,D2,DD1,DD2,Frame->True,
        Prolog->AbsolutePointSize[4],
        DisplayFunction:>$DisplayFunction,
        FrameLabel->{"P","P",Title,""},PlotRange->All];
    Return[aux1]
    ];
    
   PPPlotGEV[Sample_,lambda_,delta_,kappa_,Tail_,Title_]:=
    Module[{aux,SS=Sort[Sample],n=Length[Sample],aux1,i,alpha1},
      aux=
        Table[{i/(n+1),
            If[Tail==1,FGEVM[SS[[i]],lambda,delta,kappa],
              FGEVm[SS[[i]],lambda,delta,kappa]]},{i,1,n}];
      DD1=
        ListPlot[aux,Prolog->AbsolutePointSize[5],Axes->False,
          DisplayFunction:>Identity,PlotRange->All];
      DD2=Graphics[Line[{{0,0},{1,1}}]];
      aux1=
        Show[DD1,DD2,Frame->True,Prolog->AbsolutePointSize[4],
          DisplayFunction:>$DisplayFunction,
          FrameLabel->{"P","P",Title,""},PlotRange->All];
      Return[aux1]];

QQPlotGEV[Sample_,lambda_,delta_,kappa_,Tail_,Title_]:=
  Module[{aux,SS=Sort[Sample],n=Length[Sample],aux1},
    aux=Table[{SS[[i]],
          If[Tail==1,FGEVMInv[i/(n+1),lambda,delta,kappa],
            FGEVmInv[i/(n+1),lambda,delta,kappa]]},{i,1,n}];
    DD1=ListPlot[aux,Prolog->AbsolutePointSize[4],Axes->False,
        DisplayFunction:>Identity,PlotRange->All];
    DD2=Graphics[Line[{{SS[[1]],SS[[1]]},{SS[[n]],SS[[n]]}}]];
    aux1=Show[DD1,DD2,Frame->True,
        DisplayFunction:>$DisplayFunction,
        FrameLabel->{"Q","Q",Title,""},PlotRange->All];
    Return[aux1]
    ];
    
      
      PPPlotGPD[Sample_,alpha_,kappa_,Tail_,Title_]:=
    Module[{aux,SS=Sort[Sample],n=Length[Sample],aux1,i,alpha1},
      aux=Table[{i/(n+1),
            If[Tail==1,FGPDM[SS[[i]],alpha,kappa],
              FGPDm[SS[[i]],alpha,kappa]]},{i,1,n}];
      DD1=
        ListPlot[aux,Prolog->AbsolutePointSize[4],Axes->False,
          DisplayFunction:>Identity,PlotRange->All];
      DD2=Graphics[Line[{{0,0},{1,1}}]];
      aux1=
        Show[DD1,DD2,Frame->True,Prolog->AbsolutePointSize[4],
          DisplayFunction:>$DisplayFunction,
          FrameLabel->{"P","P",Title,""},PlotRange->All];
      Return[aux1]];
      
      
     PPPlotGPDBands[Sample_,alpha_,kappa_,Tail_,Title_]:=
  Module[{aux,SS=Sort[Sample],n=Length[Sample],aux1,i,alpha1},
    alpha1=FindAlpha[Length[Sample]];
    aux=Table[{i/(n+1),
          If[Tail==1,FGPDM[SS[[i]],alpha,kappa],
            FGPDm[SS[[i]],alpha,kappa]]},{i,1,n}];
    aux1={{0,0}};
    Do[AppendTo[
        aux1,{i/(n+1),Quantile[BetaDistribution[i,n-i+1.],alpha1/2]}],{i,1,
        n}];
    AppendTo[aux1,{1,1}];
    D1=ListPlot[aux1,
        PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
        PlotJoined->True,DisplayFunction->Identity,PlotRange->All];
    aux1={{0,0}};
    Do[AppendTo[
        aux1,{i/(n+1),Quantile[BetaDistribution[i,n-i+1.],1-alpha1/2]}],{i,1,
        n}];
    AppendTo[aux1,{1,1}];
    D2=ListPlot[aux1,
        PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
        PlotJoined->True,DisplayFunction->Identity,PlotRange->All];
    DD1=ListPlot[aux,Prolog->AbsolutePointSize[4],Axes->False,
        DisplayFunction:>Identity,PlotRange->All];
    DD2=Graphics[Line[{{0,0},{1,1}}]];
    aux1=Show[D1,D2,DD1,DD2,Frame->True,
        Prolog->AbsolutePointSize[4],
        DisplayFunction:>$DisplayFunction,
        FrameLabel->{"P","P",Title,""},PlotRange->All];
    Return[aux1]];

QQPlotGPD[Sample_,alpha_,kappa_,Tail_,threshold_,Title_]:=
    Module[{aux,SS=Sort[Sample],n=Length[Sample],aux1},
      aux=Table[{threshold+SS[[i]],
            If[Tail==1,threshold+FGPDMInv[i/(n+1),alpha,kappa],
              threshold+FGPDmInv[i/(n+1),alpha,kappa]]},{i,1,n}];
      DD1=
        ListPlot[aux,Prolog->AbsolutePointSize[4],Axes->False,
          DisplayFunction:>Identity,PlotRange->All];
      DD2=
        Graphics[
          Line[{{threshold+SS[[1]],threshold+SS[[1]]},{threshold+SS[[n]],
                threshold+SS[[n]]}}]];
      aux1=
        Show[DD1,DD2,Frame->True,
          DisplayFunction:>$DisplayFunction,
          FrameLabel->{"Q","Q",Title,""},PlotRange->All];
      Return[aux1]];

PlotGEVM[lambda_,delta_,kappa_,p1_:0.01,p2_:0.99]:=Module[{xmin,xmax}, 
    xmin=FGEVMInv[p1,lambda,delta,kappa];
    xmax=FGEVMInv[p2,lambda,delta,kappa];
    Plot[fGEVM[x,lambda,delta,kappa],{x,xmin,xmax},Frame->True,
      PlotRange->All,Axes->False]
    ];

PlotGEVm[lambda_,delta_,kappa_,p1_:0.01,p2_:0.99]:=Module[{xmin,xmax}, 
    xmin=FGEVmInv[p1,lambda,delta,kappa];
    xmax=FGEVmInv[p2,lambda,delta,kappa];
    Plot[fGEVm[x,lambda,delta,kappa],{x,xmin,xmax},Frame->True,
      PlotRange->All,Axes->False]
    ];

PlotGPDM[alpha_,kappa_,p1_:0.01,p2_:0.99]:=Module[{xmin,xmax}, 
    xmin=FGPDMInv[p1,alpha,kappa];
    xmax=FGPDMInv[p2,alpha,kappa];
    Plot[fGPDM[x,alpha,kappa],{x,xmin,xmax},Frame->True,
      PlotRange->All,Axes->False]
    ];

PlotGPDm[alpha_,kappa_,p1_:0.01,p2_:0.99]:=Module[{xmin,xmax}, 
    xmin=FGPDmInv[p1,alpha,kappa];
    xmax=FGPDmInv[p2,alpha,kappa];
    Plot[fGPDm[x,alpha,kappa],{x,xmin,xmax},Frame->True,
      PlotRange->All,Axes->False]
    ];

CompleteGEVAnalysis[Sample_,Tail_,Title_,PW_:0]:=
    
    Module[{sol,LikGumbel,LikGEV,lambdaG,deltaG,lambdaGEV,deltaGEV,kappaGEV,
        n=Length[Sample],TLR,TLRA,Z,d,TWald,i,p,q,m11,m12,m13,m22,m23,m33,M,
        M1},
      sol=EstimateGEVML[Sample,Tail];
      {LikGumbel,LikGEV}={sol[[1]],sol[[3]]};
      {lambdaG,deltaG}=sol[[2]];
      {lambdaGEV,deltaGEV,kappaGEV}=sol[[4]];
      Var=sol[[5]];
      se=sol[[6]];
      Z=Map[(#-lambdaG)/deltaG&,Sample];
      d=Sum[Z[[i]]-0.5*(Z[[i]]^2)*(1-Exp[-Z[[i]]]),{i,1,n}];
      TLR=-2*(LikGEV-LikGumbel);
      TLRA=(1-2.8/n)*TLR;
      HH1=PPPlotGEV[Sample,lambdaGEV,deltaGEV,kappaGEV,Tail,Title];
      HH2=QQPlotGEV[Sample,lambdaGEV,deltaGEV,kappaGEV,Tail,Title];
      (*Prescott and Walden asymptotic covariance matrix*)
      
      psi[u_]=
        D[Log[Gamma[u]],u];
      p=(1-kappaGEV)^2*Gamma[1-2*kappaGEV];
      q=Gamma[2-kappaGEV]*(psi[1-kappaGEV]-(1-kappaGEV)/kappaGEV);
      m11=n*p/deltaGEV^2;
      m22=n*(1-2*Gamma[2-kappaGEV]+p)/(deltaGEV^2*kappaGEV^2);
      m33=
        n*(Pi^2/6+(1-0.5772157-1/kappaGEV)^2+2*q/kappaGEV+p/kappaGEV^2)/
            kappaGEV^2;
      m12=n*(p-Gamma[2-kappaGEV])/(deltaGEV^2*kappaGEV);
      m13=-n*(q+p/kappaGEV)/(deltaGEV*kappaGEV);
      m23=
        n*(1-0.5772157-(1-Gamma[2-kappaGEV])/kappaGEV-q-
                p/kappaGEV)/(deltaGEV*
                kappaGEV^2);
      M={{m11,m12,m13},{m12,m22,m23},{m13,m23,m33}};
      M1=Inverse[M];
      TWald=kappaGEV^2/(n*M1[[3,3]]);
            If[PW!=0,
                
        Print["Prescott and Walden asymptotic covariance matrix : \r ",
                    MatrixForm[M1]]];
            Print["TLR=",TLR," p-value=",1-CDF[ChiDistribution[1],TLR]];
            Print["TLRA=",TLRA," p-value=",1-CDF[ChiDistribution[1],TLRA]];
            
      Print["TWald=",TWald," p-value=",1-CDF[ChiDistribution[1],TWald]];
            Show[GraphicsArray[{{HH1,HH2}}],GraphicsSpacing->0,
                ImageSize->{500,Automatic}];
              
      Return[{LikGumbel,LikGEV,kappaGEV,TLR,1-CDF[ChiDistribution[1],TLR],
          TLRA,1-CDF[ChiDistribution[1],TLRA],TWald,
          1-CDF[ChiDistribution[1],TWald],M1}]];
 
 Curvature[Sample_,Tail_]:=
  Module[{i,j,m,n=Length[Sample],Sample1,Sigma10,Sigma01,Sigma11,
      Sigma20,S1,S2},
      If[Tail==1,Sample1=Sort[Sample],Sample1=Sort[-Sample]];
    i=n+1-Floor[2*Sqrt[n]];j=n+1-Floor[Sqrt[n]];
    m=j-i+1;
    Sigma10=Sum[-Log[-Log[(k-0.5)/n]],{k,i,j}];
    Sigma01=Sum[Sample1[[k]],{k,i,j}];
    Sigma11=Sum[-Sample1[[k]]*Log[-Log[(k-0.5)/n]],{k,i,j}];
    Sigma20=Sum[(-Log[-Log[(k-0.5)/n]])^2,{k,i,j}];
    S1=(m*Sigma11-Sigma10*Sigma01)/(m*Sigma20-Sigma10^2);
    i=n+1-Floor[Sqrt[n]];j=n;
    m=j-i+1;
    Sigma10=Sum[-Log[-Log[(k-0.5)/n]],{k,i,j}];
    Sigma01=Sum[Sample1[[k]],{k,i,j}];
    Sigma11=Sum[-Sample1[[k]]*Log[-Log[(k-0.5)/n]],{k,i,j}];
    Sigma20=Sum[(-Log[-Log[(k-0.5)/n]])^2,{k,i,j}];
    S2=(m*Sigma11-Sigma10*Sigma01)/(m*Sigma20-Sigma10^2);
    Return[S1/S2]
    ];
    
CurvatureP[n_,NSimulations_,Tail_,cval_]:=Module[{SS={},i,sol,Sample,aux},
      Do[Sample=Table[-Log[-Log[Random[]]],{i,1,n}];
        If[Tail==0,Sample=-Sample];
        AppendTo[SS,Curvature[Sample,Tail]],{i,1,NSimulations}];
      sol=Sort[Append[SS,cval]];
      aux=(Position[sol,cval][[1,1]]-1)/NSimulations;
      Return[{aux,Min[aux,1-aux]}]];
     
fOrderStatistic[x_,r_,n_]:=f[x]*(F[x]^(r-1)*(1-F[x])^(n-r))/Beta[r,n-r+1];

FOrderStatistic[x_,r_,n_]:=r*Binomial[n,r]*Beta[F[x],r,n-r+1];

fTwoOrderStatistic[x_,y_,n_,r_,s_]:=
  If[x<y,(n!*f[x]*f[y]*
          F[x]^(r-1)*(F[y]-F[x])^(s-r-1)*(1-F[y])^(n-s))/((r-1)!*(s-r-1)!*(n-
                s)!),0];
                
PlotOrderStatistic[fx_,Fx_,r_,n_,xmin_,xmax_]:=
    Plot[(fx)*((Fx)^(r-1)*(1-(Fx))^(n-r))/Beta[r,n-r+1],{x,xmin,xmax},
      Frame->True];

PlotTwoOrderStatistics[fx_,Fx_,n_,r_,s_,xmin_,xmax_,ymin_,ymax_]:=
  Module[{f1,F1,f2,F2,HH},
    HH[x_,y_]:=((n!*(f1)*(f2)*(F1)^(r-1)*((F2)-(F1))^(s-r-1)*(1-(F2))^(n-
                      s))/((r-1)!*(s-r-1)!*(n-s)!))/.{f1->fx,F1->Fx,
          f2->fx/.{x->y},F2->Fx/.{x->y}};
    Print[HH[x,y]];
    Plot3D[If[x<y,HH[x,y],0],{x,xmin,xmax},{y,ymin,ymax},PlotRange->All,
      PlotPoints->40];
    ];

        
  RangeEstimateGPDML[Sample_,u0_,u1_,m_,Tail_]:=
    Module[{Sample1,Sample2,val={},u,aux,A1,A2,A3,B1,B2,B3,D1,D2,D3,D4,D5,D6,
        D7,D8,D9,mm,iii,val3={}},Sample2=Sort[Sample];
      Do[If[Tail==1,Sample1=Select[Sample2,#>u&]-u,
          Sample1=Select[Sample2,#<u&];Sample1=-(Sample1-u)];
        mm=Sum[Sample1[[iii]],{iii,1,Length[Sample1]}]/Length[Sample1];
        AppendTo[val3,{u,mm}];
        aux=EstimateGPDML[Sample1,1,0];
        AppendTo[val,{u,aux[[2]],Sqrt[aux[[3,1,1]]],Sqrt[aux[[3,2,2]]]}],{u,
          u0,u1,(u1-u0)/m}];
      A1=Table[{val[[i,1]],val[[i,2,1]]},{i,1,Length[val]}];
      A2=Table[{val[[i,1]],val[[i,2,1]]-val[[i,3]]},{i,1,Length[val]}];
      A3=Table[{val[[i,1]],val[[i,2,1]]+val[[i,3]]},{i,1,Length[val]}];
      B1=Table[{val[[i,1]],val[[i,2,2]]},{i,1,Length[val]}];
      B2=Table[{val[[i,1]],val[[i,2,2]]-val[[i,4]]},{i,1,Length[val]}];
      B3=Table[{val[[i,1]],val[[i,2,2]]+val[[i,4]]},{i,1,Length[val]}];
      D1=ListPlot[A1,PlotJoined->True,DisplayFunction->Identity];
      D2=ListPlot[A2,PlotJoined->True,
          PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
          DisplayFunction->Identity];
      D3=ListPlot[A3,PlotJoined->True,
          PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
          DisplayFunction->Identity];
      D4=ListPlot[B1,PlotJoined->True,DisplayFunction->Identity];
      D5=ListPlot[B2,PlotJoined->True,
          PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
          DisplayFunction->Identity];
      D6=ListPlot[B3,PlotJoined->True,
          PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
          DisplayFunction->Identity];
      D7=Show[D1,D2,D3,Frame->True,PlotRange->All,
          FrameLabel->{"u","lambda"},Axes->False,
          DisplayFunction->$DisplayFunction];
      D8=Show[D4,D5,D6,Frame->True,PlotRange->All,
          FrameLabel->{"u","kappa"},Axes->False,
          DisplayFunction->$DisplayFunction];
      D9=ListPlot[val3,PlotJoined->True,DisplayFunction->Identity,
          Frame->True,FrameLabel->{"u","E[X-u|X>u]"}];
      Show[GraphicsArray[{{D7},{D8},{D9}},ImageSize->300]]];
      
     
   
      
ToUnitFrechetMarginal[Sample_,Tail_,options___]:=
    Module[{optionslist,Estim,Gumb,sol,lambdaGEV,deltaGEV,kappaGEV,Sample1},
      optionslist={options};
      {Estim,Gumb}={Estimates,Gumbel}/.optionslist/.Options[\
ToUnitFrechetMarginal];
      If[Tail==0,Sample1=-Sample,Sample1=Sample];
      If[Gumb,sol=EstimateGumbelML[Sample,1,0];{lambdaGEV,deltaGEV,kappaGEV}={\
sol[[2]],sol[[3]],0},sol=EstimateGEVML[Sample,1,0];{lambdaGEV,deltaGEV,\
kappaGEV}=sol[[4]]];
      If[Gumb,Sample1=Map[Exp[(#-lambdaGEV)/deltaGEV]&,Sample1],
        Sample1=Map[(1-kappaGEV*(#-lambdaGEV)/deltaGEV)^(-1/kappaGEV)&,
            Sample1]];
      If[Tail==0,Sample1=-Sample1;lambdaGEV=-lambdaGEV];
      If[Estim,Return[{Sample1,{lambdaGEV,deltaGEV,kappaGEV}}],
        Return[Sample1]]];
     
FLogistic[x_,y_,alpha_]:=Exp[-(x^(-1/alpha)+y^(-1/alpha))^alpha];

fLogistic[x_,y_,alpha_]=Simplify[D[FLogistic[x,y,alpha],x,y]];

LKLogistic[Sample_,
    alpha_]:=-Sum[
      Log[fLogistic[Sample[[i,1]],Sample[[i,2]],alpha]],{i,1,Length[Sample]}];\


EstimateLogisticML[Sample_,p_:1]:=Module[{sol,alpha0},
    sol=FindMinimum[LKLogistic[Sample,alpha],{alpha,0.00001,1}];
    alpha0=alpha/.sol[[2]];
    If[p==1,Print["Likelihood=",sol[[1]]," alpha=",alpha0]];
    Return[{sol[[1]],alpha0}]
    ];
    
EmpiricalTwodimF[Sample_,x_,
      y_]:=(Sum[
            If[((Sample[[i,1]]<=x ) && (Sample[[i,2]]<=y)),1.0,0.0],{i,1,
              Length[Sample]}]+0.5)/(Length[Sample]+1);
  
EmpiricalTwodimF[Sample_,x_,
      y_]:=(Sum[
            If[((Sample[[i,1]]<=x ) && (Sample[[i,2]]<=y)),1.0,0.0],{i,1,
              Length[Sample]}]+0.5)/(Length[Sample]+1);

EstimateLogisticLS[Sample_,p_:1]:=
  Module[{pp,n=Length[Sample],QQ,i,alpha0,sol},
    pp=Table[EmpiricalTwodimF[Sample,Sample[[i,1]],Sample[[i,2]]],{i,1,n}];
    QQ[alpha_]:=
      Sum[(FLogistic[Sample[[i,1]],Sample[[i,2]],alpha]-
                pp[[i]])^2/(pp[[i]]*(1-pp[[i]])/Length[Sample]),{i,1,
          Length[Sample]}];
    sol=FindMinimum[QQ[alpha],{alpha,0.3,0.4}];
    alpha0=alpha/.sol[[2]];
    If[p==1,Print["SSQ=",-sol[[1]]," alpha=",alpha0]];
    Return[{sol[[1]],alpha0}]
    ];

PointProcMLLogistic[Sample_,p_:1]:=Module[{i,W,sol,alpha0,PPMLE},
      W=Table[
          Sample[[i,1]]/(Sample[[i,1]]+Sample[[i,2]]),{i,1,Length[Sample]}];
      PPMLE[
          alpha_]:=-Sum[
            Log[1/alpha-1]-(1+1/alpha)*Log[W[[i]]*(1-W[[i]])]+(alpha-2)*
                Log[W[[i]]^(-1/alpha)+(1-W[[i]])^(-1/alpha)],{i,1,
              Length[Sample]}];
      sol=FindMinimum[PPMLE[alpha],{alpha,0.2,0.3},MaxIterations->100];
      alpha0=alpha/.sol[[2]];
      If[p==1,Print["LogLikelihood=",-sol[[1]]," alpha=",alpha0]];
      Return[{sol[[1]],alpha0}]];
               
PointProcMLDirichlet[Sample_,p_:1]:=
    Module[{i,W,sol,alpha0,beta0,PPMLE},
      W=Table[Sample[[i,1]]/(Sample[[i,1]]+Sample[[i,2]]),{i,1,
            Length[Sample]}];
      PPMLE[alpha_,
          beta_]:=-Sum[
            Log[alpha^2*beta^2*
                Gamma[alpha^2+
                    beta^2+1]*(alpha^2*
                      W[[i]])^(alpha^2-1)*(beta^2*(1-W[[i]]))^(beta^2-1)/(2*
                      Gamma[alpha^2]*
                      Gamma[beta^2]*(alpha^2*W[[i]]+
                            beta^2*(1-W[[i]]))^(alpha^2+beta^2+1))],{i,1,
              Length[Sample]}];
      sol=FindMinimum[PPMLE[alpha,beta],{alpha,0.00001,1},{beta,0.00001,5}];
      {alpha0,beta0}={alpha^2,beta^2}/.sol[[2]];
      If[p==1,Print["LogLikelihood=",-sol[[1]]," alpha=",alpha0," \
beta=",beta0],
        " beta=",beta0];
      Return[{sol[[1]],alpha0,beta0}]];
                  
      FGumbelLogistic[x_,y_,lambda1_,delta1_,lambda2_,delta2_,alpha_]:=
        Exp[-(Exp[(lambda1-x)/(delta1*alpha)]+
                  Exp[(lambda2-y)/(delta2*alpha)])^alpha];
                  
      fGumbelLogistic[x_,y_,lambda1_,delta1_,lambda2_,delta2_,alpha_]=
        D[FGumbelLogistic[x,y,lambda1,delta1,lambda2,delta2,alpha],x,y];
          
PointProcJointMLLogisticGumbel[Sample_,p_:1]:=
    Module[{lambda10,delta10,lambda20,delta20,alpha0,i,sol,sol1,sol2,Samp1,
        Samp2,MLGumbelLogistic},
      MLGumbelLogistic[lambda1_,delta1_,lambda2_,delta2_,
          alpha_]:=-Sum[
            Log[fGumbelLogistic[Sample[[i,1]],Sample[[i,2]],lambda1,delta1^2,\
lambda2,
                delta2^2,alpha^2]],{i,1,Length[Sample]}];
      Samp1=Table[Sample[[i,1]],{i,1,Length[Sample]}];
      Samp2=Table[Sample[[i,2]],{i,1,Length[Sample]}];
      sol1=EstimateGEVML[Samp1,1,0];
      sol2=EstimateGEVML[Samp2,1,0];
      {lambda10,delta10,lambda20,delta20,alpha0}={sol1[[4,1]],sol1[[4,2]],
          sol2[[4,1]],sol2[[4,2]],0.3};
      Print[
        "Initial estimate = ",{lambda10,delta10,lambda20,delta20,alpha0}];
        {delta10,delta20,alpha0}={delta10^(1/2),delta20^(1/2),alpha0^(1/2)};
      sol=
        FindMinimum[
          MLGumbelLogistic[lambda1,delta1,lambda2,delta2,alpha],{lambda1,
            lambda10,lambda10*1.1},{delta1,delta10,delta10*1.1},{lambda2,\
lambda20,lambda20*1.1},{delta2,
            delta20,delta20*1.1},{alpha,0.2,0.3},MaxIterations->100];
      {lambda10,delta10,lambda20,delta20,
          alpha0}={lambda1,delta1,lambda2,delta2,alpha}/.sol[[2]];
          {delta10,delta20,alpha0}={delta10^2,delta20^2,alpha0^2};
          If[p==1,
        Print["LogLikelihood=",-sol[[1]]," lambda1=",lambda10," delta1=",
          delta10," lambda2=",lambda20," delta2=",delta20," alpha=",alpha0]];
      Return[{sol[[1]],lambda10,delta10,lambda20,delta20,alpha0}]];
      
PlotLogistic[Sample_,lambda11_,delta11_,lambda21_,delta21_,alpha1_,Size_:400]:\
=
  Module[{D1Data,D1Plot,xmin,xmax,ymin,ymax,Samp1,Samp2,n=Length[Sample],i},
    Samp1=Table[Sample[[i,1]],{i,1,n}];
    Samp2=Table[Sample[[i,2]],{i,1,n}];
    xmin=Min[Samp1];
    xmax=Max[Samp1];
    ymin=Min[Samp2];
    ymax=Max[Samp2];
    D1Plot=
      ContourPlot[
        fGumbelLogistic[x,y,lambda11,delta11,lambda21,delta21,alpha1],{x,xmin,\
xmax},{y,
          ymin,ymax},ContourShading->False,PlotPoints->100,
        Frame->True,PlotRange->All,DisplayFunction->Identity];
    D1Data=
      ListPlot[Sample,Frame->True,Axes->False,Frame->True,
        DisplayFunction->Identity];
    Return[
      Show[D1Plot,D1Data,Frame->True,AspectRatio->1,Prolog->AbsolutePointSize[\
3],DisplayFunction:>$DisplayFunction,ImageSize->Size]];
    ];
    
    FGumbelDirichlet[x_,y_,lambda1_,delta1_,lambda2_,delta2_,alpha_,beta_]:=
        H[Exp[(x-lambda1)/delta1],Exp[(y-lambda2)/delta2],alpha,beta];
      fGumbelDirichlet[x_,y_,lambda1_,delta1_,lambda2_,delta2_,alpha_,beta_]=
        D[FGumbelDirichlet[x,y,lambda1,delta1,lambda2,delta2,alpha,beta],x,y];\

        
        
   PlotDirichlet[Sample_,lambda11_,delta11_,lambda21_,delta21_,alpha1_,beta1_]\
:=
    Module[{D1Data,D1Plot,xmin,xmax,ymin,ymax,Samp1,Samp2,n=Length[Sample],i},\

      H[x_,y_,alpha_,beta_]:=
        Exp[-(1-Beta[alpha*x/(alpha*x+beta*y),alpha+1,beta])/
              x-(1-Beta[alpha*x/(alpha*x+beta*y),alpha,beta+1])/y];
      Samp1=Table[Sample[[i,1]],{i,1,n}];
      Samp2=Table[Sample[[i,2]],{i,1,n}];
      xmin=Min[Samp1];
      xmax=Max[Samp1];
      ymin=Min[Samp2];
      ymax=Max[Samp2];
      D1Plot=
        ContourPlot[
          fGumbelDirichlet[x,y,lambda11,delta11,lambda21,delta21,alpha1,beta1]\
,{x,xmin,
            xmax},{y,ymin,ymax},ContourShading->False,
          PlotPoints->100,Frame->True,
          DisplayFunction->Identity];
      D1Data=
        ListPlot[Sample,Frame->True,Axes->False,Frame->True,
          DisplayFunction->Identity];
      Return[
        Show[D1Plot,D1Data,Frame->True,
          DisplayFunction:>$DisplayFunction]];];
  
 F[Sample1_, x_] := 
    Length[Select[Sample1, #1<= x &]]/Length[Sample1];
    
PlotHistogramGEV[Sample_,lambda_,delta_,kappa_,d_: 8]:=
    Module[{n=Length[Sample],Sample1,PSample,xmin,xmax,DS,DD,increment,aux},
      Print["d=",d];Sample1=Sort[Sample]; xmin=Min[Sample1]; 
      xmax=Max[Sample1]; increment=(xmax-xmin)/d;
      DS=Plot[fGEVM[x,lambda,delta,kappa],{x,xmin,xmax},
          PlotStyle->{AbsoluteThickness[3]},
          DisplayFunction->Identity]; DD={}; 
      Do[x1=x+increment; y=(F[Sample1,x1]-F[Sample1,x])/increment;
        AppendTo[DD,Line[{{x,0},{x,y},{x1,y},{x1,0}}]],{x,xmin,
          xmax-increment/100,increment}]; 
      aux=Show[DS,Graphics[DD],PlotRange->All,Frame->True,
          AxesOrigin->{xmin,0},Axes->True,
          FrameLabel->{"X","p"},DisplayFunction->$DisplayFunction]; 
      Return[aux]];
          
Bisection[xa_,xb_,Dj_,Errx_,Erry_]:=Module[{F0,F1,F2,x0=xa,x1=xb,x2},
      F0=h[xa,Dj];
      F1=h[xb,Dj];
      If[F0*F1<=0,While[Abs[x0-x1]>Errx,x2=(x0+x1)*0.5;F2=h[x2,Dj];
          If[F2*F0>0.0,x0=x2;F0=F2,x1=x2;F1=F2]],
        If[Abs[F0]<Erry,x2=xa,
          If[Abs[F1]<Erry,x2=xb,
            Print["No Solution F0=",F0," F1=",F1," x=",x]]]];Return[x2]];

EstimateGEVEPM1[Sample_,Tail_,Option_]:=
    Module[{Sample1,i,j,n=Length[Sample],A1n,Ajn,Aj1,C1,Cj,Cn,Errx=0.00001,
        Erry=0.00001,xa,xb,kappa,tol=0.000001,lambdaval,deltaval,lowval,nn,n1,\
n2},
      h[kappa_,Dj_]:=(1-Ajn^kappa)/(1-A1n^kappa)-Dj;
      If[Tail==1,Sample1=Sort[Sample],Sample1=Sort[-Sample]];
      C1=-Log[(1-0.35)/n];
      Cn=-Log[(n-0.35)/n];
      A1n=C1/Cn;
      kappaval=Table[0,{n-2}];
      deltaval=kappaval;
      lowval=kappaval;
      Do[Cj=-Log[(j-0.35)/n];
        Ajn=Cj/Cn;
        Aj1=Cj/C1;
        Dj=(Sample1[[j]]-Sample1[[n]])/(Sample1[[1]]-Sample1[[n]]);
        If[Dj<Log[Ajn]/Log[A1n],xa=tol;xb=Log[Dj]/Log[Aj1];
          kappa=Bisection[xa,xb,Dj,Errx,Erry],xa=Log[1-Dj]/Log[Ajn];xb=-tol;
          kappa=Bisection[xa,xb,Dj,Errx,Erry]];
        If[Abs[kappa]<0.01,delta=(Sample1[[1]]-Sample1[[j]])/(Log[Cj/C1]);
          lambda=Sample1[[1]]+delta*Log[C1],
          delta=kappa*(Sample1[[1]]-Sample1[[n]])/(Cn^kappa-C1^kappa);
          low=Sample1[[1]]+delta*C1^kappa/kappa];
        lowval[[j-1]]=low;
        kappaval[[j-1]]=kappa;
        deltaval[[j-1]]=delta
        ,{j,2,n-1}];
        nn=n-2;
      kappaval=Sort[kappaval];
      deltaval=Sort[deltaval];
      lowval=Sort[lowval];
      n1=Max[Floor[0.2*n+0.4999],1];
      n2=Min[Floor[0.8*n+0.4999],n-2];
      If[Option==1,kappa=Sum[kappaval[[j]],{j,n1,n2}]/(n2-n1+1);
      delta=Sum[deltaval[[j]],{j,n1,n2}]/(n2-n1+1);
      low=Sum[lowval[[j]],{j,n1,n2}]/(n2-n1+1),
      kappa=(kappaval[[Floor[(nn+1)/2]]]+kappaval[[Floor[(nn+2)/2]]])/2;
      delta=(deltaval[[Floor[(nn+1)/2]]]+deltaval[[Floor[(nn+2)/2]]])/2;
      low=(lowval[[Floor[(nn+1)/2]]]+lowval[[Floor[(nn+2)/2]]])/2];
      lambda=low-delta/kappa;
      If[Tail==0,lambda=-lambda];
      Return[{lambda,delta,kappa}]];
      
      EstimateGEVEPM[Sample_,Tail_,Option_]:=
  Module[{Sample1,i,j,n=Length[Sample],A1n,Ajn,Aj1,C1,Cj,Cn,Errx=0.00001,
      Erry=0.00001,xa,xb,kappa,tol=0.000001,lambdaval,deltaval,lowval,nn,n1,
      n2,nestimates=0,jjj,iii,aaa},
    h[kappa_,Dj_]:=(1-Ajn^kappa)/(1-A1n^kappa)-Dj;
    If[Tail==1,Sample1=Sort[Sample],Sample1=Sort[-Sample]];
    aaa=n*(n-1)*(n-2)/6;
    kappaval=Table[0,{aaa}];
      deltaval=kappaval;
      lowval=kappaval;
      lambdaval=kappaval;
    Do[C1=-Log[(iii-0.35)/n];
      Cn=-Log[(jjj-0.35)/n];
      A1n=C1/Cn;
      Do[If[(Sample1[[iii]]!=Sample1[[j]]) &&(Sample1[[j]]!=
                Sample1[[jjj]]), Cj=-Log[(j-0.35)/n];
          Ajn=Cj/Cn;
          Aj1=Cj/C1;
          Dj=(Sample1[[j]]-Sample1[[jjj]])/(Sample1[[iii]]-Sample1[[jjj]]);
          If[Dj<Log[Ajn]/Log[A1n],xa=tol;xb=Log[Dj]/Log[Aj1];
            kappa=Bisection[xa,xb,Dj,Errx,Erry],xa=Log[1-Dj]/Log[Ajn];
            xb=-tol;
            kappa=Bisection[xa,xb,Dj,Errx,Erry]];
          
          If[Abs[kappa]<0.01,
            delta=(Sample1[[iii]]-Sample1[[j]])/(Log[Cj/C1]);
            lambda=Sample1[[iii]]+delta*Log[C1],
            delta=kappa*(Sample1[[iii]]-Sample1[[jjj]])/(Cn^kappa-C1^kappa);
            low=Sample1[[iii]]+delta*C1^kappa/kappa];
          
          If[(kappa>0 && low>=Sample1[[iii]]) || (kappa<0 && 
                  low<= Sample1[[jjj]]),nestimates=nestimates+1];
          lambda=low-delta/kappa;
          (*Print["nestimates=",nestimates," iii=",iii," j=",j," jjj=",jjj,
              " lambda=",lambda," delta=",delta," kappa=",kappa ];*)
         
           lambdaval[[nestimates]]=lambda;
          deltaval[[nestimates]]=delta;
          kappaval[[nestimates]]=kappa;
          lowval[[nestimates]]=low]
        ,{j,iii+1,jjj-1}],{iii,1,n-2},{jjj,iii+2,n}];
    nn=nestimates;
    lambdaval=Take[lambdaval,nn];
    deltaval=Take[deltaval,nn];
    kappaval=Take[kappaval,nn];
    lowval=Take[lowval,nn];
    lambdaval=Sort[lambdaval];
    deltaval=Sort[deltaval];
    kappaval=Sort[kappaval];
    lowval=Sort[lowval];
    (*Print["lambdaval=",lambdaval];*)
    (*Print["deltaval=",
        deltaval];*)
    (*Print["kappaval=",kappaval];*)
    (*Print[
        "lowval=",lowval];*)
    n1=Max[Floor[0.2*nestimates+0.4999],1];
    n2=Min[Floor[0.8*nestimates+0.4999],nn-2];
    If[Option==1,lambda=Sum[lambdaval[[j]],{j,n1,n2}]/(n2-n1+1);
      delta=Sum[deltaval[[j]],{j,n1,n2}]/(n2-n1+1);
      kappa=Sum[kappaval[[j]],{j,n1,n2}]/(n2-n1+1);
      low=Sum[lowval[[j]],{j,n1,n2}]/(n2-n1+1),
      lambda=(lambdaval[[Floor[(nn+1)/2]]]+lambdaval[[Floor[(nn+2)/2]]])/2;
      delta=(deltaval[[Floor[(nn+1)/2]]]+deltaval[[Floor[(nn+2)/2]]])/2;
      kappa=(kappaval[[Floor[(nn+1)/2]]]+kappaval[[Floor[(nn+2)/2]]])/2;
      low=(lowval[[Floor[(nn+1)/2]]]+lowval[[Floor[(nn+2)/2]]])/2];
    (*lambda=low-delta/kappa;*)
    If[Tail==0,lambda=-lambda];
    Print["nestimates=",nestimates," low=",low," delta=",delta," kappa=",
      kappa," lambda=",lambda];
    Return[{lambda,delta,kappa}]];
      
    EstimateGEVEPMThreshold[TruncatedSample_,Tail_,n_,Option_]:=
    Module[{Sample1,i,j,n0=Length[TruncatedSample],A1n,Ajn,Aj1,C1,Cj,Cn,Errx=\
0.0001,
        Erry=0.0001,xa,xb,kappa,tol=0.000001,lambdaval,deltaval,lowval,nn,n1,\
n2},
      h[kappa_,Dj_]:=(1-Ajn^kappa)/(1-A1n^kappa)-Dj;
      If[Tail==1,Sample1=Sort[TruncatedSample],Sample1=Sort[-TruncatedSample]]\
;
      C1=-Log[(n-n0+1-0.35)/n];
      Cn=-Log[(n-0.35)/n];
      A1n=C1/Cn;
      nn=n0-2;
      kappaval=Table[0,{nn}];
      deltaval=Table[0,{nn}];
      lowval=Table[0,{nn}];
      Do[
        Cj=-Log[(n-n0+j-0.35)/n];
        Ajn=Cj/Cn;
        Aj1=Cj/C1;
        Dj=(Sample1[[j]]-Sample1[[n0]])/(Sample1[[1]]-Sample1[[n0]]);
        If[Dj<Log[Ajn]/Log[A1n],xa=tol;xb=Log[Dj]/Log[Aj1];
          kappa=Bisection[xa,xb,Dj,Errx,Erry],xa=Log[1-Dj]/Log[Ajn];xb=-tol;
          kappa=Bisection[xa,xb,Dj,Errx,Erry]];
        If[Abs[kappa]<0.01,delta=(Sample1[[1]]-Sample1[[j]])/(Log[Cj/C1]);
          lambda=Sample1[[1]]+delta*Log[C1],
          delta=kappa*(Sample1[[1]]-Sample1[[n0]])/(Cn^kappa-C1^kappa);
          low=Sample1[[1]]+C1^kappa/kappa];
        lowval[[j-1]]=low;
        kappaval[[j-1]]=kappa;
        deltaval[[j-1]]=delta;,{j,2,n0-1}];
      kappaval=Sort[kappaval];
      deltaval=Sort[deltaval];
      lowval=Sort[lowval];
      n1=Max[Floor[0.2*n0+0.4999],1];
      n2=Min[Floor[0.8*n0+0.4999],n0-2];
      If[Option==1,kappa=Sum[kappaval[[j]],{j,n1,n2}]/(n2-n1+1);
      delta=Sum[deltaval[[j]],{j,n1,n2}]/(n2-n1+1);
      low=Sum[lowval[[j]],{j,n1,n2}]/(n2-n1+1),
      kappa=(kappaval[[Floor[(nn+1)/2]]]+kappaval[[Floor[(nn+2)/2]]])/2;
      delta=(deltaval[[Floor[(nn+1)/2]]]+deltaval[[Floor[(nn+2)/2]]])/2;
      low=(lowval[[Floor[(nn+1)/2]]]+lowval[[Floor[(nn+2)/2]]])/2];
      lambda=low-delta/kappa;
      If[Tail==0,lambda=-lambda];
      Return[{lambda,delta,kappa}]];
      
PlotPDFOrderStatistic[Fx_,var_,r_,n_,xmin_,xmax_,p_:1]:=
    Module[{f,F,aux,lplot={},r1=Flatten[List[r]],i,rr},
      F=Function[var,Fx];
      f[x_]=D[F[x],x];
      Do[rr=r1[[i]];
        aux=Plot[
            f[x]*(F[x]^(rr-1)*(1-F[x])^(n-rr))/Beta[rr,n-rr+1],{x,xmin,xmax},
            Frame->True,PlotRange->All,
            DisplayFunction->Identity];
        AppendTo[lplot,aux],{i,1,Length[r1]}];
      If[p==1,
        aux=Show[#,DisplayFunction->$DisplayFunction]&@@{lplot},aux=Show[#,\
DisplayFunction->Identity]&@@{lplot}];
      Return[aux]];
      
PlotPDFOrderStatistic1[fx_,Fx_,var_,r_,n_,xmin_,xmax_,p_:1]:=
    Module[{f,F,aux,lplot={},r1=Flatten[List[r]],i,rr},f=Function[var,fx];
      F=Function[var,Fx];
      Do[rr=r1[[i]];
        aux=Plot[
            f[x]*(F[x]^(rr-1)*(1-F[x])^(n-rr))/Beta[rr,n-rr+1],{x,xmin,xmax},
            Frame->True,PlotRange->All,
            DisplayFunction->Identity];
        AppendTo[lplot,aux],{i,1,Length[r1]}];
      If[p==1,
        aux=Show[#,DisplayFunction->$DisplayFunction]&@@{lplot},aux=Show[#,\
DisplayFunction->Identity]&@@{lplot}];
      Return[aux]];

PlotCDFOrderStatistic[Fx_,var_,r_,n_,xmin_,xmax_,p_:1]:=
    Module[{F,aux,lplot={},r1=Flatten[List[r]],i,rr},F=Function[var,Fx];
      Do[rr=r1[[i]];
        aux=Plot[rr*Binomial[n,rr]*Beta[F[x],rr,n-rr+1],{x,xmin,xmax},
            Frame->True,PlotRange->All,
            DisplayFunction->Identity];
        AppendTo[lplot,aux],{i,1,Length[r1]}];
      If[p==1,
        aux=Show[#,DisplayFunction->$DisplayFunction]&@@{lplot},aux=Show[#,\
DisplayFunction->Identity]&@@{lplot}];
      Return[aux]];

PlotPDFTwoOrderStatistics[Fx_,var_,r_,s_,n_,xmin_,xmax_,ymin_,ymax_,p_:1]:=
    Module[{f,F1,f2,F2,HH,aux},F=Function[var,Fx];
    f[x_]=D[F[x],x];
      HH[x_,y_]:=
        If[x<y,((n!*f[x]*f[y]*
                  F[x]^(r-1)*(F[y]-F[x])^(s-r-1)*(1-F[y])^(n-s))/((r-1)!*(s-
                        r-1)!*(n-s)!)),0];
      If[p==1,
        aux=Plot3D[HH[x,y],{x,xmin,xmax},{y,ymin,ymax},PlotRange->All,
            PlotPoints->40,DisplayFunction->$DisplayFunction],
        aux=
          Plot3D[HH[x,y],{x,xmin,xmax},{y,ymin,ymax},PlotRange->All,
            PlotPoints->40,DisplayFunction->Identity]
        ];Return[aux]];
        
PlotPDFTwoOrderStatistics1[fx_,Fx_,var_,r_,s_,n_,xmin_,xmax_,ymin_,ymax_,p_:1]\
:=
    Module[{f1,F1,f2,F2,HH,aux},F=Function[var,Fx];f=Function[var,fx];
      HH[x_,y_]:=
        If[x<y,((n!*f[x]*f[y]*
                  F[x]^(r-1)*(F[y]-F[x])^(s-r-1)*(1-F[y])^(n-s))/((r-1)!*(s-
                        r-1)!*(n-s)!)),0];
      If[p==1,
        aux=Plot3D[HH[x,y],{x,xmin,xmax},{y,ymin,ymax},PlotRange->All,
            PlotPoints->40,DisplayFunction->$DisplayFunction],
        aux=
          Plot3D[HH[x,y],{x,xmin,xmax},{y,ymin,ymax},PlotRange->All,
            PlotPoints->40,DisplayFunction->Identity]
        ];Return[aux]];
          
EstimateGPDCH1[Sample_,Tail_,Option_,p_:0]:=
    Module[{Sample1,i,j,n=Length[Sample],Errx=0.00001,Erry=0.00001,xa,xb,\
kappa,
        alpha,alphak,tol=0.000001,lambdaval,lowval,h,n1,n2},
      Bisection1[xa_,xb_,Errx_,Erry_]:=
        Module[{F0,F1,F2,x0=xa,x1=xb,x2},F0=h[xa];
          F1=h[xb];
          If[F0*F1<=0,
            While[Abs[x0-x1]>Errx,x2=(x0+x1)*0.5;F2=h[x2];
              If[F2*F0>0.0,x0=x2;F0=F2,x1=x2;F1=F2]],
            If[Abs[F0]<Erry,x2=xa,
              If[Abs[F1]<Erry,x2=xb,
                Print["No Solution F0=",F0," F1=",F1," x=",x]]]];Return[x2]];
      h[kappa_]:=
        Sample1[[i]]*(1-(1-j/(n+1))^kappa)-
          Sample1[[j]]*(1-(1-i/(n+1))^kappa);
      If[Tail==1,Sample1=Sort[Sample],Sample1=Sort[-Sample]];
      kappaval=Table[0,{n-1}];
      lowval=kappaval;
      count=0;
      Do[Do[count+=1;
          d0=Sample1[[i]]/Sample1[[j]]-Log[1-i/(n+1)]/Log[1-j/(n+1)];
          If[d0>0,xa=tol;xb=Log[1-Sample1[[i]]/Sample1[[j]]]/Log[1-i/(n+1)];
            kappa=Bisection1[xa,xb,Errx,Erry],
            xa=Log[Sample1[[i]]/Sample1[[j]]]/Log[(1-i/(n+1))/(1-j/(n+1))];
            xb=-tol;
            kappa=Bisection1[xa,xb,Errx,Erry]];
          alphak=Sample1[[j]]/(1-(1-j/(n+1))^kappa);
          lowval[[count]]=alphak;
          kappaval[[count]]=kappa;     
          If[p==1,
            Print["kappa=",kappa," alpha/k=",alphak," alpha=",
              alphak*kappa]],{j,n,n}],{i,1,n-1}];
      kappaval=Sort[kappaval];
      lowval=Sort[lowval];
      n1=Max[Floor[0.2*n+0.4999],1];
      n2=Min[Floor[0.8*n+0.4999],n-1];
      If[Option==1,kappa=Sum[kappaval[[j]],{j,n1,n2}]/(n2-n1+1);
      low=Sum[lowval[[j]],{j,n1,n2}]/(n2-n1+1),
      kappa=(kappaval[[Floor[(count+1)/2]]]+kappaval[[Floor[(count+2)/2]]])/2;\

      low=(lowval[[Floor[(count+1)/2]]]+lowval[[Floor[(count+2)/2]]])/2];
      alpha=Abs[low*kappa];
      Return[{alpha,kappa}]];
      
    
    
  
    
EstimateGPDCH[Sample_,Tail_,Option_,p_:0]:=
  Module[{Sample1,i,j,n=Length[Sample],Errx=0.00001,Erry=0.00001,xa,xb,kappa,
      alpha,alphak,tol=0.000001,lambdaval,lowval,h,n1,n2,aaa,count,maxVal,
      error=-1000},
    Bisection1[xa_,xb_,Errx_,Erry_]:=
      Module[{F0,F1,F2,x0=xa,x1=xb,x2=xa},F0=h[xa];
        F1=h[xb];
        If[F0*F1<=0,While[Abs[x0-x1]>Errx,x2=(x0+x1)*0.5;F2=h[x2];
            If[F2*F0>0.0,x0=x2;F0=F2,x1=x2;F1=F2]],
          If[Abs[F0]<Erry,x2=xa,
            If[Abs[F1]<Erry,x2=xb,
              Print["No Solution F0=",F0," F1=",F1," x=",x]]]];Return[x2]];
    h[kappa_]:=
      Sample1[[i]]*(1-(1-j/(n+1))^kappa)-Sample1[[j]]*(1-(1-i/(n+1))^kappa);
    If[Tail==1,Sample1=Sort[Sample],Sample1=Sort[-Sample]];
    maxVal=Sample1[[n]];
    aaa=n*(n-1)/2;
    kappaval=Table[0,{aaa}];
    lowval=kappaval;
    count=0;
    Do[Do[
        If[Sample1[[i]]!=Sample1[[j]],
          d0=Sample1[[i]]/Sample1[[j]]-Log[1-i/(n+1)]/Log[1-j/(n+1)];
          
          If[d0>0,xa=tol;
            xb=Log[1-Sample1[[i]]/Sample1[[j]]]/Log[1.0-i/(n+1)],
            xa=Log[Sample1[[i]]/Sample1[[j]]]/Log[(1.0-i/(n+1))/(1-j/(n+1))];
            xb=-tol];
          kappa=Bisection1[xa,xb,Errx,Erry];
          error=Max[error,Abs[h[kappa]]];
          alphak=kappa*Sample1[[j]]/(1-(1-j/(n+1))^kappa);
          If[alphak>0,count=count+1;lowval[[count]]=alphak;
            kappaval[[count]]=kappa;     
            
            If[p==1,
              Print["kappa=",kappa," alpha/k=",alphak," alpha=",
                alphak*kappa]]]]
        ,{j,i+1,n}],{i,1,n-1}];
    kappaval=Take[kappaval,count];
    lowval=Take[lowval,count];
    kappaval=Sort[kappaval];
    lowval=Sort[lowval];
    n1=Max[Floor[0.2*count+0.4999],1];
    n2=Min[Floor[0.8*count+0.4999],count-1];
    Print["count=",count," n1=",n1," n2=",n2];
    If[Option==1,kappa=Sum[kappaval[[j]],{j,n1,n2}]/(n2-n1+1);
      low=Sum[lowval[[j]],{j,n1,n2}]/(n2-n1+1),
      kappa=(kappaval[[Floor[(count+1)/2]]]+kappaval[[Floor[(count+2)/2]]])/
          2;
      low=(lowval[[Floor[(count+1)/2]]]+lowval[[Floor[(count+2)/2]]])/2];
    alpha=low;
    Return[{alpha,kappa}];
    ]
    
 EstimateGPDLS[Sample_,Tail_]:=
  Module[{Sample1,i,n=Length[Sample],lambda0,kappa0,lambda,kappa},
    h[lambda_,kappa_]:=
      If[Abs[kappa]>0.01,
        Sum[(Sample1[[i]]-lambda*(1-(1-(i-0.5)/n)^kappa)/kappa)^2,{i,1,n}],
        Sum[(Sample1[[i]]+lambda*Log[1-(i-0.5)/n])^2,{i,1,n}]];
    If[Tail==1,Sample1=Sort[Sample],Sample1=Sort[-Sample]];
    sol=EstimateGPDPWM[Sample1,1];
    {lambda0,kappa0}={sol[[1]],sol[[2]]};
    sol=FindMinimum[
        h[lambda,kappa],{lambda,lambda0,lambda0*1.1},{kappa,kappa0,
          kappa0*1.1}];
    Return[{lambda,kappa}/.sol[[2]]];]
    
 EstimateGPDPWM[Sample_,Tail_]:=
    Module[{Sample1,t,mean,n=Length[Sample],Var},
      If[Tail==0,Sample1=Sort[-Sample],Sample1=Sort[Sample]];
      mean=Sum[Sample1[[i]],{i,1,n}]/n;
      t=Sum[(1-(i-0.35)/n)*Sample1[[i]],{i,1,n}]/n;
      kappa=(4*t-mean)/(mean-2*t);
      alpha=2*mean*t/(mean-2*t);
      Var=(1/(n*(1+2*kappa)*(3+2*kappa)))*{{kappa^2*(7+18*kappa+11*kappa^2+2*
                      kappa^3),
              kappa*(2+kappa)*(2+6*kappa+7*kappa^2+2*kappa^3)},{kappa*(2+
                    kappa)*(2+6*kappa+7*kappa^2+2*kappa^3),(1+
                    kappa)*(2+kappa)^2*(1+kappa+2*kappa^2)}};
      Return[{alpha,kappa,Var,{Sqrt[Var[[1,1]]],Sqrt[Var[[2,2]]]}}]
      ];
     
MakeBivariateSample[MarginalX_,MarginalY_]:=Transpose[{MarginalX,MarginalY}];
     

PlotBivariateSample[Sample_,Size_:600]:=
  ListPlot[Sample,Frame->True,Axes->False,ImageSize->Size,AspectRatio->1,\
Prolog->AbsolutePointSize[4],PlotRange->All];
   
   
DomainofAttraction[xp_,var_]:=Module[{Q,aux,kappamax,kappamin},
    Q=Function[var,xp];
    aux=Limit[(Q[epsilon]-Q[2*epsilon])/(Q[2*epsilon]-Q[4*epsilon]),
        epsilon->0];kappamin=-Log[aux]/Log[2];
    aux=Limit[(Q[1-epsilon]-Q[1-2*epsilon])/(Q[1-2*epsilon]-Q[1-4*epsilon]),
        epsilon->0];kappamax=-Log[aux]/Log[2];
    Return[N[{kappamin,kappamax}]]
    ];
    
EstimateGEVLS[Sample_,Tail_]:=
  Module[{Sample1,i,n=Length[Sample],lambda0,delta0,kappa0,lambda,delta,
      kappa},h[lambda_,delta_,kappa_]:=
      Sum[(Sample1[[i]]-lambda-delta*(1-(-Log[i/(n+1)])^kappa)/kappa)^2,{i,1,
          n}];
    If[Tail==1,Sample1=Sort[Sample],Sample1=Sort[-Sample]];
    sol=EstimateGEVPWM[Sample1,1];
    {lambda0,delta0,kappa0}=sol;
    sol=FindMinimum[
        h[lambda,delta,kappa],{lambda,lambda0,lambda0*1.1},{delta,delta0,
          delta0*1.1},{kappa,kappa0,kappa0*1.1}];
    {lambda,delta,kappa}={lambda,delta,kappa}/.sol[[2]];
    If[Tail==0,lambda=-lambda];
    Return[{lambda,delta,kappa}];]
    
EstimateGEVPWM[Sample_,Tail_]:=
    Module[{Sample1,b0,b1,b2,Dj,xa,xb,Errx=0.0001,Erry=0.0001,
        n=Length[Sample]},
        If[Tail==0,Sample1=-Sample,Sample1=Sample];
      Sample1=Sort[Sample1];
      b0=Sum[Sample1[[i]],{i,1,n}]/n;
      b1=Sum[((i-0.35)/n)*Sample1[[i]],{i,1,n}]/n;
      b2=Sum[((i-0.35)/n)^2*Sample1[[i]],{i,1,n}]/n;
      Dj=(3*b2-b0)/(2*b1-b0);
      If[Dj<1,Print["No solution for the PWM method"];Return[]];
      h[kappa_,Dj_]:=(1-3^(-kappa))/(1-2^(-kappa))-Dj;
      xa=-7;xb=10;
      kappa=Bisection[xa,xb,Dj,Errx,Erry];
      delta=kappa*(2*b1-b0)/(Gamma[1+kappa]*(1-2^(-kappa)));
      lambda=b0-delta*(1-Gamma[1+kappa])/kappa;
      If[Tail==0,lambda=-lambda];
      Return[{lambda,delta,kappa}]];
                  
     FGEVLogistic[x_,y_,lambda1_,delta1_,kappa1_,lambda2_,delta2_,kappa2_,\
alpha_]:=
        Exp[-(If[Abs[kappa1]<0.002,Exp[(lambda1-x)/(delta1*alpha)],(1-kappa1*(\
x-lambda1)/delta1)^(1/(alpha*kappa1))]
        +If[Abs[kappa2]<0.002,Exp[(lambda2-y)/(delta2*alpha)],(1-
                        kappa2*(y-lambda2)/delta2)^(1/(alpha*kappa2))])^
              alpha];
          
      fGEVLogistic[x_,y_,lambda1_,delta1_,kappa1_,lambda2_,delta2_,kappa2_,\
alpha_]=
        Simplify[
          D[FGEVLogistic[x,y,lambda1,delta1,kappa1,lambda2,delta2,kappa2,\
alpha],x,
            y]];
            
   PointProcJointMLLogisticGEV[Sample_,p_:1]:=
    Module[{lambda10,delta10,kappa10,lambda20,delta20,kappa20,alpha0,i,sol,
        sol1,sol2,Samp1,Samp2,MLGEVLogistic,minval1,maxval1,minval2,
        maxval2},
      MLGEVLogistic[lambda1_,delta1_,kappa1_,lambda2_,delta2_,kappa2_,alpha_]:\
=
        If[(((kappa1>0) && (lambda1+delta1/kappa1)>
                      maxval1)|| ((kappa1<0) && (lambda1+delta1/kappa1)<
                      minval1)) && (((kappa2>0) && (lambda2+delta2/kappa2)>
                      maxval2)|| ((kappa2<0) && (lambda2+delta2/kappa2)<
                      minval2)),-Sum[
              Log[fGEVLogistic[Sample[[i,1]],Sample[[i,2]],lambda1,delta1,\
kappa1,
                  lambda2,delta2,kappa2,alpha]],{i,1,Length[Sample]}],10^6];
      Samp1=Table[Sample[[i,1]],{i,1,Length[Sample]}];
      Samp2=Table[Sample[[i,2]],{i,1,Length[Sample]}];
      minval1=Min[Samp1];
      maxval1=Max[Samp1];
      minval2=Min[Samp2];
      maxval2=Max[Samp2];
      sol1=EstimateGEVML[Samp1,1,0];
      sol2=EstimateGEVML[Samp2,1,0];
      {lambda10,delta10,kappa10,lambda20,delta20,kappa20,alpha0}={sol1[[4,1]],\

          sol1[[4,2]],sol1[[4,3]],sol2[[4,1]],sol2[[4,2]],sol2[[4,3]],0.3};
      Print[
        "Initial estimate = ",{lambda10,delta10,kappa10,lambda20,delta20,
          kappa20,alpha0}];
      sol=
        FindMinimum[
          MLGEVLogistic[lambda1,delta1,kappa1,lambda2,delta2,kappa2,alpha],{\
lambda1,
            lambda10,lambda10*1.1},{delta1,delta10,delta10*1.1},{kappa1,
            kappa10,kappa10*1.1},{lambda2,lambda20,lambda20*1.1},{delta2,
            delta20,delta20*1.1},{kappa2,kappa20,kappa20*1.1},{alpha,0.2,0.3},\

          MaxIterations->100];
      {lambda10,delta10,kappa10,lambda20,delta20,kappa20,
          alpha0}={lambda1,delta1,kappa1,lambda2,delta2,kappa2,
            alpha}/.sol[[2]];
      If[p==1,
        Print["LogLikelihood=",-sol[[1]]," lambda1=",lambda10," delta1=",
          delta10," kappa1=",kappa10," lambda2=",lambda20," delta2=",delta20,
          " kappa2=",kappa20," alpha=",alpha0]];
      Return[{sol[[1]],lambda10,delta10,kappa10,lambda20,delta20,kappa20,
          alpha0}]];
          
PlotGEVLogistic[Sample_,lambda1_,delta1_,kappa1_,lambda2_,delta2_,kappa2_,
          alpha_,Size_:400]:=Module[{D1,D2,aux,minval,maxval},
    minval=Map[Min[#]& , Transpose[Sample]];
    maxval=Map[Max[#]& , Transpose[Sample]];
    D1=ListPlot[Sample,Frame->True,Axes->False,PlotRange->All,
        DisplayFunction->Identity];
    D2=ContourPlot[
        fGEVLogistic[x,y,lambda1,delta1,kappa1,lambda2,delta2,kappa2,
          alpha],{x,minval[[1]],
          maxval[[1]]},{y,minval[[2]],maxval[[2]]},ContourShading->False,
        PlotRange->All,PlotPoints->100,
        DisplayFunction->Identity];
    aux=Show[D1,D2,DisplayFunction->$DisplayFunction,AspectRatio->1,Prolog->\
AbsolutePointSize[3],ImageSize->Size];
    Return[aux]
    ];

Beta1[u_,alpha_,beta_]:=
    Gamma[alpha+beta]*Beta[u,alpha,beta]/Gamma[alpha]/Gamma[beta];

FDirichlet[x_,y_,alpha_,beta_]:=
    Exp[-(1-Beta1[alpha*x/(alpha*x+beta*y),alpha+1,beta])/
          x-(Beta1[alpha*x/(alpha*x+beta*y),alpha,beta+1])/y];

fDirichlet[x_,y_,alpha_,beta_]=If[x>0 && y>0,f1Dirichlet[x,y,alpha,beta],0];

f1Dirichlet[x_,y_,alpha_,beta_]=Simplify[D[FDirichlet[x,y,alpha,beta],x,y]];

PlotGEVDirichlet[Sample_,lambda1_,delta1_,kappa1_,lambda2_,delta2_,kappa2_,
      alpha_,beta_,Size_:400]:=
    Module[{D1,D2,aux,minval,maxval},minval=Map[Min[#]&,Transpose[Sample]];
      maxval=Map[Max[#]&,Transpose[Sample]];
      D1=ListPlot[Sample,Frame->True,Axes->False,
          PlotRange->All,DisplayFunction->Identity];
      D2=ContourPlot[
          fDirichlet[(1-kappa1*(x-lambda1)/delta1)^(1/(alpha*kappa1)),(1-
                  kappa2*(y-lambda2)/delta2)^(1/(alpha*kappa2)),alpha,
            beta],{x,minval[[1]],maxval[[1]]},{y,minval[[2]],maxval[[2]]},
          ContourShading->False,PlotRange->All,PlotPoints->100,
          DisplayFunction->Identity];
      aux=
        Show[D1,D2,DisplayFunction->$DisplayFunction,AspectRatio->1,
          Prolog->AbsolutePointSize[4],ImageSize->Size];
      Return[aux]];


EstimateGEVDirichletML[Sample_,p_:1]:=
    Module[{lambda10,delta10,kappa10,lambda20,delta20,kappa20,alpha0,beta0,i,
        sol,sol1,sol2,Samp1,Samp2,aux,minval1,maxval1,minval2,maxval2},
        MLGEVDirichlet[lambda1_,delta1_,kappa1_,lambda2_,delta2_,kappa2_,\
alpha_,
      beta_]:=Module[{alpha1,beta1,aux},
      If[((((kappa1>0)&&(lambda1+delta1/kappa1)>
                      maxval1)||((kappa1<0)&&(lambda1+delta1/kappa1)<
                      minval1))&&(((kappa2>0)&&(lambda2+delta2/kappa2)>
                      maxval2)||((kappa2<0)&&(lambda2+delta2/kappa2)<
                      minval2)) && alpha>0 && beta>0),
        If[alpha>100,alpha1=100-Exp[-(alpha-100)]+1,alpha1=alpha];
        If[beta>100,beta1=100-Exp[-(beta-100)]+1,beta1=beta];
        aux=-Sum[
              Log[fDirichlet[(1-
                        kappa1*(Sample[[i,1]]-lambda1)/delta1)^(1/(alpha1*
                            kappa1)),(1-
                        kappa2*(Sample[[i,2]]-lambda2)/delta2)^(1/(alpha1*
                            kappa2)),alpha1,beta1]],{i,1,Length[Sample]}],
        aux=10^6];
    (*  Print["aux>=",aux,
        " Estimates: ",{lambda1,delta1,kappa1,lambda2,delta2,kappa2,alpha,
          beta}];*)
      Return[aux]];
      Samp1=Table[Sample[[i,1]],{i,1,Length[Sample]}];
      Samp2=Table[Sample[[i,2]],{i,1,Length[Sample]}];
      minval1=Min[Samp1];
      maxval1=Max[Samp1];
      minval2=Min[Samp2];
      maxval2=Max[Samp2];
      sol1=EstimateGEVML[Samp1,1,0];
      sol2=EstimateGEVML[Samp2,1,0];
      {lambda10,delta10,kappa10,lambda20,delta20,kappa20,alpha0,
          beta0}={sol1[[4,1]],sol1[[4,2]],sol1[[4,3]],sol2[[4,1]],sol2[[4,2]],\

          sol2[[4,3]],100,100};
      Print[
        "Initial estimate = ",{lambda10,delta10,kappa10,lambda20,delta20,
          kappa20,alpha0,beta0}];
      sol=
        FindMinimum[
          MLGEVDirichlet[lambda1,delta1,kappa1,lambda2,delta2,kappa2,alpha,
            beta],{lambda1,lambda10,lambda10*1.1},{delta1,delta10,
            delta10*1.1},{kappa1,kappa10,kappa10*1.1},{lambda2,lambda20,
            lambda20*1.1},{delta2,delta20,delta20*1.1},{kappa2,kappa20,
            kappa20*1.1},{alpha,alpha0,alpha0*1.1},{beta,beta0,beta0*1.1},
          MaxIterations->100];
      {lambda10,delta10,kappa10,lambda20,delta20,kappa20,alpha0,
          beta0}={lambda1,delta1,kappa1,lambda2,delta2,kappa2,alpha,
            beta}/.sol[[2]];
      If[p==1,
        Print["LogLikelihood=",-sol[[1]]," lambda1=",lambda10," delta1=",
          delta10," kappa1=",kappa10," lambda2=",lambda20," delta2=",delta20,
          " kappa2=",kappa20," alpha=",alpha0," beta=",beta0]];
      Return[{sol[[1]],lambda10,delta10,kappa10,lambda20,delta20,kappa20,
          alpha0,beta0}]];
          
Redond[x_]:=
    Module[{n=Floor[Log[x]/Log[10.]]},Return[{Floor[x/(10^n)]*(10^n),n}]];

Ticks1[xmin_,xmax_]:=Module[{amin,amax,tt,k},amin=Redond[xmin];
      amax=Redond[xmax];
      k=1;
      While[Length[tt]<6,tt={};
        Do[AppendTo[tt,j*10.0^i],{i,amin[[2]],amax[[2]]},{j,1,10,k}];k=k/2;
        tt=Union[tt];
        tt=Select[tt,#>=xmin&&#<=xmax&]];
      Return[N[tt]]];

Ticks2[xmin_,xmax_]:=Module[{amin,amax,tt={},k,ord1={1,5,2,8,4,6,3,9,7}},
      amin=Redond[xmin];
      amax=Redond[xmax];
      If[amin[[2]]==amax[[2]],Return[Ticks1[xmin,xmax]],
        Do[
          AppendTo[tt,ord1[[j]]*10.0^i],{j,1,Length[ord1]},{i,amin[[2]],
            amax[[2]]}];
        tt=Select[tt,#>=xmin&&#<=xmax&];
        Return[N[tt]]]];
        
Ticks3[xmin_,xmax_]:=
    Module[{aux},
      If[xmin>0&&xmax>0,aux=Ticks2[xmin,xmax],
        If[xmin<0&&xmax<0,aux=-Ticks2[-xmax,-xmin],
          If[xmin<0&&xmax>0,
            aux=Join[-Ticks2[10^(Redond[-xmin][[2]]-1),-xmin],
                Ticks2[10^(Redond[xmax][[2]]-1),xmax]];PrependTo[aux,0]]]];
      Return[aux]];

TicksMarks[xmin_,xmax_]:=Module[{amin,amax,tt,n1,n2},
      amin=Redond[xmin];
      amax=Redond[xmax];
        tt={};
        Do[AppendTo[tt,j*10.0^i],{i,amin[[2]],amax [[2]]},{j,0,10,5}];
      tt=Union[tt];
      tt=Select[tt,#>=xmin && #<=xmax &];
      Return[tt]
      ];
      
      
ProbabilityPaperLine2[Sample_,paper_,lambda_,NpointsLine_,options___]:=
    Module[{XTran,YTran,PlotPointPosition,
        Probabilities={0.00001,0.0001,0.001,0.005,0.01,0.02,0.05,0.1,0.2,0.3,
            0.4,0.5,0.6,0.7,0.8,0.9,0.95,0.98,0.99,0.995,0.999,0.9999,
            0.99999},Probabilities0,Probabilities1,ReturnPeriods,i,i1,i2,iii,
        val,a,b,NpointsLine1=NpointsLine,SampleSize=Length[Sample],ProbPoints,\

        Points,Sample1,Y,meanX,meanY,DD1,DD2,DD3,DD8,DD9,SXX,SXY,Error=0,
        Tlog={2,3,5,6,8},XTransf={1,2,3,1,4,4,1,3},
        label={"X","X","X","X","X","X","X","X"},MX,bb,Bandplus,Bandminus,aux5,\

        aux6,u,sigma,alpha=5/100,x1,x2,y1,y2,xmin,xmax,ymin,ymax,auxx,Smin,
        Smax,Grid1,Grid2,tt,a0,a1,delta,ttfin,valid,optionslist},
           optionslist={options};
      PlotSize=PlotSize/.optionslist/.Options[ProbabilityPaperLine2];
      If[(Intersection[Distr,{paper}]=={} && 
            Intersection[{paper},Range[Length[Distr]]]=={}),
        Print["Unknown paper"];Error=1;Return["Error"]];
      If[StringQ[paper],paper1=Position[Distr,paper][[1,1]],paper1=paper];
      
            If[paper1==3||paper1==8,
        If[lambda<Max[Sample],
          Print["****** Erroneous lambda. Minimum value of lambda must be: ",
            Max[Sample]," ******"];Error=1]];
      If[paper1==5||paper1==6,
        If[Min[Sample]<lambda,
          Print["****** Erroneous lambda. Maximum value of lambda must be: ",
            Min[Sample]," ******"];Error=1;]];
            
      i1=-1;
      Do[If[SampleSize/(SampleSize+1)<Probabilities[[i]],i1=i1+1],{i,1,
          Length[Probabilities]}];
      Probabilities=Take[Probabilities,{i1+1,Length[Probabilities]-i1}];
      PTran[k_,x_]:=
        If[k<3,Quantile[NormalDistribution[0,1],x],
          If[k<6,-Log[-Log[x]],If[k<9,Log[-Log[1-x]],0]]];
      PTran1[k_,x_]:=
        If[k<3,CDF[NormalDistribution[0,1],x],
          If[k<6,Exp[-Exp[-x]],If[k<9,1-Exp[-Exp[x]],0]]];
      XTran[k_,x_]:=
        Which[k==1,x,k==2,x,k==3,Log[1/(lambda-x)],
          k==4,x,k==5,Log[x-lambda],k==6,Log[x-lambda],
          k==7,x,k==8,Log[1/(lambda-x)],True,0];
      Smin=Min[Sample];
      Smax=Max[Sample];
      tt=Ticks3[Smin,Smax];
      ttfin={};
	a0=XTran[paper1,Smin];
	a1=XTran[paper1,Smax];
	delta=(a1-a0)/20;
	Do[valid=True;
	  Do[If[Abs[XTran[paper1,tt[[i]]]-XTran[paper1,ttfin[[j]]]]<delta,valid=\
False]
 	   ,{j,1,Length[ttfin]}];
	  If[valid,AppendTo[ttfin,tt[[i]]]]
 	 ,{i,1,Length[tt]}];
	ttfin=Sort[ttfin];
      Grid1=Map[{N[XTran[paper1,#]],ToString[#]}&,ttfin];
      Grid2=Map[N[XTran[paper1,#]]&,ttfin];
      PlotPointPosition[i_,SampleSize_]:=N[i/(SampleSize+1)];
      If[NpointsLine1==0,NpointsLine1=SampleSize];
      If[paper1==3||paper1==8,
        label[[paper1]]=
          StringJoin["X  Transformation: -log(",ToString[lambda],
            "- X)"]];
      If[lambda>0,
        If[paper1==5||paper1==6,
          label[[paper1]]=
            StringJoin["X  Transformation: log(X - ",ToString[lambda],
              ")"]]];
      If[lambda<0,
        If[paper1==5||paper1==6,
          label[[paper1]]=
            StringJoin["X  Transformation: log(X ",ToString[lambda],
              ")"]]];
      Print["Plotting a ",Distr[[paper1]]," probability paper"];
      Sample1=Sort[Sample];

      If[Abs[NpointsLine1]>SampleSize||Abs[NpointsLine1]<2,
        Print["Error in the number of points to be used for the regression \
line it must be 2 <= n <= ",SampleSize];Error=1];
      If[Error!=1,Sample1=Map[XTran[paper1,#]&,Sample1];
        If[NpointsLine1>0,i1=SampleSize-NpointsLine1+1;i2=SampleSize,i1=1;
          i2=-NpointsLine1];
        NpointsLine1=Abs[NpointsLine1];
        Probabilities0=Probabilities;
        Probabilities1=Map[{N[PTran[paper1,#]],ToString[#]}&,Probabilities0];
        ReturnPeriods=
          If[paper1<6,
            Table[{Probabilities1[[i,1]],
                Floor[100/(1-Probabilities0[[i]])+0.5]/100.0},{i,1,
                Length[Probabilities1]}],
            Table[{Probabilities1[[i,1]],
                Floor[100/Probabilities0[[i]]+0.5]/100.0},{i,1,
                Length[Probabilities1]}]];
        ProbPoints=
          N[Map[PTran[paper1,#]&,
              Table[PlotPointPosition[i,SampleSize],{i,1,SampleSize}]]];
        Points=Table[{Sample1[[i]],ProbPoints[[i]]},{i,1,SampleSize}];
        (*Regression line*)
        (*Print["Sample1=",Sample1];*)
        
        meanX=Sum[Sample1[[i]],{i,i1,i2}]/NpointsLine1;
        meanY=Sum[ProbPoints[[i]],{i,i1,i2}]/NpointsLine1;
        SXX=Sum[(Sample1[[i]]-meanX)^2,{i,i1,i2}]/NpointsLine1;
        SXY=
          Sum[(Sample1[[i]]-meanX)*(ProbPoints[[i]]-meanY),{i,i1,i2}]/
            NpointsLine1;
        b=SXY/SXX;a=meanY-b*meanX;
        If[a>0,Print["Regression line: v=",b,"u+",a],Print["Regression line: \
v=",b,"u",a]];
        (*Print[{{Sample1[[1]],a+b*Sample1[[1]]},{Sample1[[SampleSize]],
                  a+b*Sample1[[SampleSize]]}}];*)(*Print["i1=",i1," i2=",i2,
              " NpointsLine=",NpointsLine1];*)
        
        ExtremesLine={{Sample1[[1]],a+b*Sample1[[1]]},{Sample1[[SampleSize]],
              a+b*Sample1[[SampleSize]]}};
        ymax=Probabilities1[[Length[Probabilities1],1]];
        ymin=Probabilities1[[1,1]];
        If[ymax<ExtremesLine[[2,2]],auxx=(ymax-a)/b;
          ExtremesLine[[2]]={auxx,ymax}];
        If[ymin>ExtremesLine[[1,2]],auxx=(ymin-a)/b;
          ExtremesLine[[1]]={auxx,ymin}];
        MX=Table[{1,Sample1[[i]]},{i,i1,i2}];
        bb=Inverse[Transpose[MX].MX];
        sigma=
          Sqrt[Sum[(ProbPoints[[i]]-a-b*Sample1[[i]])^2,{i,i1,i2}]/(i2-i1)];
        (*Print["sigma=",sigma," bb=",bb];*)
        Bandplus={};
        Bandminus={};
        xmin=Sample1[[1]];
        xmax=Sample1[[SampleSize]];
        ymin=Probabilities1[[1,1]];
        ymax=Probabilities1[[Length[Probabilities1],1]];
        Do[
          iii=(SampleSize+1)*PTran1[paper1,a+b*u];
          
          If[iii<SampleSize+0.95&&iii>0.05,
            aux5=PTran[paper1,
                Quantile[BetaDistribution[iii,SampleSize-iii+1],1-alpha/2]];
            If[u>=xmin&&u<=xmax&&
                aux5>=ymin&&aux5<=ymax,
              AppendTo[Bandplus,{u,aux5}]];
            
            aux6=PTran[paper1,
                Quantile[BetaDistribution[iii,SampleSize-iii+1],alpha/2]];
            
            If[u>=xmin&&u<=xmax&&
                aux6>=ymin&&aux6<=ymax,
              AppendTo[Bandminus,{u,aux6}]]],{u,Sample1[[i1]],
            Sample1[[i2]],(Sample1[[i2]]-Sample1[[i1]])/100}];
        (*Print["a=",a," b=",b];
          Print["val=",val];
          Print["Points=",Points];Print["Probabilities0=",Probabilities0];
          Print["ProbPoints=",ProbPoints];
          Print["ReturnPeriods=",ReturnPeriods];
          Print["Probabilities1=",Probabilities1];
          Print["ExtremesLine=",ExtremesLine];
          Print["Range=",{{Points[[1,1]],
                Points[[1,2]]},{Points[[SampleSize,1]],
                Points[[SampleSize,2]]}}];*)
        
        If[Intersection[Tlog,{paper1}]!={},
          DD8=ListPlot[Bandplus,
              FrameTicks->{Grid1,Probabilities1,Automatic,ReturnPeriods},
              Frame->True,
              GridLines->{Grid2,PTran[paper1,Probabilities]},
              Prolog->AbsolutePointSize[4],
              FrameLabel->{label[[paper1]],"Probability",Distr[[paper1]],
                  "Return Period"},PlotJoined->True,
              PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
              DisplayFunction:>Identity];
          
          DD9=ListPlot[Bandminus,
              FrameTicks->{Grid1,Probabilities1,Automatic,ReturnPeriods},
              Frame->True,
              GridLines->{Grid2,PTran[paper1,Probabilities]},
              Prolog->AbsolutePointSize[4],
              FrameLabel->{label[[paper1]],"Probability",Distr[[paper1]],
                  "Return Period"},PlotJoined->True,
              PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
              DisplayFunction:>Identity];
          
          DD1=ListPlot[ExtremesLine,
              FrameTicks->{Grid1,Probabilities1,Automatic,ReturnPeriods},
              Frame->True,
              GridLines->{Grid2,PTran[paper1,Probabilities]},
              Prolog->AbsolutePointSize[4],
              FrameLabel->{label[[paper1]],"Probability",Distr[[paper1]],
                  "Return Period"},PlotJoined->True,
              PlotStyle->AbsoluteThickness[1],
              DisplayFunction:>Identity];
          
          DD2=ListPlot[Points,
              FrameTicks->{Grid1,Probabilities1,Automatic,ReturnPeriods},
              GridLines->{Grid2,PTran[paper1,Probabilities]},
              Prolog->AbsolutePointSize[4],
              FrameLabel->{label[[paper1]],"Probability",Distr[[paper1]],
                  "Return Period"},DisplayFunction:>Identity],
          DD8=ListPlot[Bandplus,
              FrameTicks->{Automatic,Probabilities1,Automatic,
                  ReturnPeriods},
              AxesOrigin->{Points[[1,1]],Points[[1,2]]},Frame->True,
              GridLines->{Automatic,PTran[paper1,Probabilities]},
              Prolog->AbsolutePointSize[4],
              FrameLabel->{"X","Probability",Distr[[paper1]],
                  "Return Period"},PlotJoined->True,
              PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
              DisplayFunction:>Identity];
          
          DD9=ListPlot[Bandminus,
              FrameTicks->{Automatic,Probabilities1,Automatic,
                  ReturnPeriods},
              AxesOrigin->{Points[[1,1]],Points[[1,2]]},Frame->True,
              GridLines->{Automatic,PTran[paper1,Probabilities]},
              Prolog->AbsolutePointSize[4],
              FrameLabel->{"X","Probability",Distr[[paper1]],
                  "Return Period"},PlotJoined->True,
              PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
              DisplayFunction:>Identity];
          
          DD1=ListPlot[ExtremesLine,
              FrameTicks->{Automatic,Probabilities1,Automatic,
                  ReturnPeriods},
              AxesOrigin->{Points[[1,1]],Points[[1,2]]},Frame->True,
              GridLines->{Automatic,PTran[paper1,Probabilities]},
              Prolog->AbsolutePointSize[4],
              FrameLabel->{"X","Probability",Distr[[paper1]],
                  "Return Period"},PlotJoined->True,
              PlotStyle->AbsoluteThickness[1],
              DisplayFunction:>Identity];
          
          DD2=ListPlot[Points,
              PlotRange->{{Points[[1,1]],
                    Points[[1,2]]},{Points[[SampleSize,1]],
                    Points[[SampleSize,2]]}},
              FrameTicks->{Automatic,Probabilities1,Automatic,
                  ReturnPeriods},
              AxesOrigin->{Points[[1,1]],Points[[1,2]]},
              GridLines->{Automatic,PTran[paper1,Probabilities]},
              Prolog->AbsolutePointSize[4],
              FrameLabel->{"X","Probability",Distr[[paper1]],
                  "Return Period"},DisplayFunction:>Identity]],
        Return[]];
      DD3=
        Show[DD1,DD2,DD8,DD9,Axes->False,Frame->True,
          DisplayFunction:>$DisplayFunction,PlotRange->All,
          ImageSize->{PlotSize,Automatic}];
      Return[DD3]]; 
      
ProbabilityPaperLine[Sample_,paper_,lambda_,NpointsLine_,options___]:=
    Module[{XTran,YTran,PlotPointPosition,Probabilities,
        Probabilities11={0.00001,0.0001,0.001,0.005,0.01,0.02,0.05,0.1,0.2,
            0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.95,0.98,0.99,0.995,0.999,0.9999,
            0.99999},Probabilities0,Probabilities1,ReturnPeriods,i,i1,i2,iii,
        val,a,b,NpointsLine1=NpointsLine,SampleSize=Length[Sample],ProbPoints,\

        Points,Sample1,Y,meanX,meanY,DD1,DD2,DD3,DD8,DD9,SXX,SXY,Error=0,
        Tlog={2,3,5,6,8},XTransf={1,2,3,1,4,4,1,3},
        label={"X","X","X","X","X","X","X","X"},MX,bb,Bandplus,Bandminus,aux5,\

        aux6,u,sigma,alpha=5/100,x1,x2,y1,y2,xmin,xmax,ymin,ymax,auxx,Smin,
        Smax,Grid1,Grid2,tt,a0,a1,delta,ttfin,valid,Gr1,optionslist,contr=0,\
imSize},
      optionslist={options};
      Gr1=Grids/.optionslist/.Options[ProbabilityPaperLine];
      Probabilities=
        ProbabilityTicks/.optionslist/.Options[ProbabilityPaperLine];
      If[!ListQ[Probabilities],Probabilities=Probabilities11];
      imSize=PlotSize/.optionslist/.Options[ProbabilityPaperLine];
      If[(Intersection[Distr,{paper}]=={}&&
            Intersection[{paper},Range[Length[Distr]]]=={}),
        Print["Unknown paper"];Error=1;Return["Error"]];
      If[StringQ[paper],paper1=Position[Distr,paper][[1,1]],paper1=paper];
      If[paper1==3||paper1==8,
        If[lambda<Max[Sample],
          Print["****** Erroneous lambda. Minimum value of lambda must be: ",
            Max[Sample]," ******"];Error=1]];
      If[paper1==5||paper1==6,
        If[Min[Sample]<lambda,
          Print["****** Erroneous lambda. Maximum value of lambda must be: ",
            Min[Sample]," ******"];Error=1;]];
      i1=-1;
      Do[If[SampleSize/(SampleSize+1)<Probabilities[[i]],i1=i1+1],{i,1,
          Length[Probabilities]}];
      Probabilities=Take[Probabilities,{i1+1,Length[Probabilities]-i1}];
      PTran[k_,x_]:=
        If[k<3,Quantile[NormalDistribution[0,1],x],
          If[k<6,-Log[-Log[x]],If[k<9,Log[-Log[1-x]],0]]];
      PTran1[k_,x_]:=
        If[k<3,CDF[NormalDistribution[0,1],x],
          If[k<6,Exp[-Exp[-x]],If[k<9,1-Exp[-Exp[x]],0]]];
      XTran[k_,x_]:=
        Which[k==1,x,k==2,x,k==3,Log[1/(lambda-x)],
          k==4,x,k==5,Log[x-lambda],k==6,Log[x-lambda],
          k==7,x,k==8,Log[1/(lambda-x)],True,0];       
      Smin=Min[Sample];
      Smax=Max[Sample];
      tt=Gr1;
      If[((!ListQ[Gr1])&&(Intersection[Tlog,{paper1}]!={})),
        tt=Ticks3[Smin,Smax];contr=1];
      If[((!ListQ[Gr1])&&(Intersection[Tlog,{paper1}]=={})),
        Grid1=Automatic;Grid2=Automatic,
        If[!ListQ[Gr1],tt=tt,
          tt=Select[Gr1,#>=Smin&&#<=Smax&]];
        ttfin={};
        a0=XTran[paper1,Smin];
        a1=XTran[paper1,Smax];
        delta=(a1-a0)/20;
        If[contr==1,Do[valid=True;
            
            Do[If[Abs[XTran[paper1,tt[[i]]]-XTran[paper1,ttfin[[j]]]]<delta,
                valid=False],{j,1,Length[ttfin]}];
            If[valid,AppendTo[ttfin,tt[[i]]]],{i,1,Length[tt]}],ttfin=tt];
        ttfin=Sort[ttfin];
        Grid1=Map[{N[XTran[paper1,#]],ToString[#]}&,ttfin];
        Grid2=Map[N[XTran[paper1,#]]&,ttfin];];
      PlotPointPosition[i_,SampleSize_]:=N[i/(SampleSize+1)];
      If[NpointsLine1==0,NpointsLine1=SampleSize];
      If[paper1==3||paper1==8,
        label[[paper1]]=
          StringJoin["X  Transformation: - log(",ToString[lambda]," - X)"]];
      If[lambda>0,
        If[paper1==5||paper1==6,
          label[[paper1]]=
            StringJoin["X  Transformation: log(X - ",ToString[lambda],")"]]];
      If[lambda<0,
        If[paper1==5||paper1==6,
          label[[paper1]]=
            StringJoin["X  Transformation: log(X ",ToString[lambda],")"]]];
          Print["Plotting a ",Distr[[paper1]]," probability paper"];
      Sample1=Sort[Sample];
      If[Abs[NpointsLine1]>SampleSize||Abs[NpointsLine1]<2,
        Print["Error in the number of points to be used for the regression \
line it must be 2 <= n <= ",SampleSize];Error=1];
      If[Error!=1,Sample1=Map[XTran[paper1,#]&,Sample1];
        If[NpointsLine1>0,i1=SampleSize-NpointsLine1+1;i2=SampleSize,i1=1;
          i2=-NpointsLine1];
        NpointsLine1=Abs[NpointsLine1];
        Probabilities0=Probabilities;
        Probabilities1=Map[{N[PTran[paper1,#]],ToString[#]}&,Probabilities0];
        ReturnPeriods=
          If[paper1<6,
            Table[{Probabilities1[[i,1]],
                Floor[100/(1-Probabilities0[[i]])+0.5]/100.0},{i,1,
                Length[Probabilities1]}],
            Table[{Probabilities1[[i,1]],
                Floor[100/Probabilities0[[i]]+0.5]/100.0},{i,1,
                Length[Probabilities1]}]];
        ProbPoints=
          N[Map[PTran[paper1,#]&,
              Table[PlotPointPosition[i,SampleSize],{i,1,SampleSize}]]];
        Points=Table[{Sample1[[i]],ProbPoints[[i]]},{i,1,SampleSize}];
        (*Regression line*)(*Print["Sample1=",Sample1];*)meanX=
          Sum[Sample1[[i]],{i,i1,i2}]/NpointsLine1;
        meanY=Sum[ProbPoints[[i]],{i,i1,i2}]/NpointsLine1;
        SXX=Sum[(Sample1[[i]]-meanX)^2,{i,i1,i2}]/NpointsLine1;
        SXY=
          Sum[(Sample1[[i]]-meanX)*(ProbPoints[[i]]-meanY),{i,i1,i2}]/
            NpointsLine1;
        b=SXY/SXX;a=meanY-b*meanX;
        If[a>0,Print["Regression line: v=",b,"u+",a],Print["Regression line: \
v=",b,"u",a]];
        (*Print[{{Sample1[[1]],a+b*Sample1[[1]]},{Sample1[[SampleSize]],
                  a+b*Sample1[[SampleSize]]}}];*)(*Print["i1=",i1," i2=",i2,
              " NpointsLine=",NpointsLine1];*)
          ExtremesLine={{Sample1[[1]],
              a+b*Sample1[[1]]},{Sample1[[SampleSize]],
              a+b*Sample1[[SampleSize]]}};
        ymax=Probabilities1[[Length[Probabilities1],1]];
        ymin=Probabilities1[[1,1]];
        If[ymax<ExtremesLine[[2,2]],auxx=(ymax-a)/b;
          ExtremesLine[[2]]={auxx,ymax}];
        If[ymin>ExtremesLine[[1,2]],auxx=(ymin-a)/b;
          ExtremesLine[[1]]={auxx,ymin}];
        MX=Table[{1,Sample1[[i]]},{i,i1,i2}];
        bb=Inverse[Transpose[MX].MX];
        sigma=
          Sqrt[Sum[(ProbPoints[[i]]-a-b*Sample1[[i]])^2,{i,i1,i2}]/(i2-i1)];
        (*Print["sigma=",sigma," bb=",bb];*)Bandplus={};
        Bandminus={};
        xmin=Sample1[[1]];
        xmax=Sample1[[SampleSize]];
        ymin=Probabilities1[[1,1]];
        ymax=Probabilities1[[Length[Probabilities1],1]];
        Do[iii=(SampleSize+1)*PTran1[paper1,a+b*u];
          
          If[iii<SampleSize+0.95&&iii>0.05,
            aux5=PTran[paper1,
                Quantile[BetaDistribution[iii,SampleSize-iii+1],1-alpha/2]];
            
            If[u>=xmin&&u<=xmax&&
                aux5>=ymin&&aux5<=ymax,
              AppendTo[Bandplus,{u,aux5}]];
            
            aux6=PTran[paper1,
                Quantile[BetaDistribution[iii,SampleSize-iii+1],alpha/2]];
            
            If[u>=xmin&&u<=xmax&&
                aux6>=ymin&&aux6<=ymax,
              AppendTo[Bandminus,{u,aux6}]]],{u,Sample1[[i1]],
            Sample1[[i2]],(Sample1[[i2]]-Sample1[[i1]])/100}];
        (*Print["a=",a," b=",b];
          Print["val=",val];
          Print["Points=",Points];Print["Probabilities0=",Probabilities0];
          Print["ProbPoints=",ProbPoints];
          Print["ReturnPeriods=",ReturnPeriods];
          Print["Probabilities1=",Probabilities1];
          Print["ExtremesLine=",ExtremesLine];
          
          Print["Range=",{{Points[[1,1]],
                Points[[1,2]]},{Points[[SampleSize,1]],
                Points[[SampleSize,2]]}}];*)If[
          Intersection[Tlog,{paper1}]!={},
          DD8=ListPlot[Bandplus,
              FrameTicks->{Grid1,Probabilities1,Automatic,ReturnPeriods},
              Frame->True,
              GridLines->{Grid2,PTran[paper1,Probabilities]},
              Prolog->AbsolutePointSize[4],
              FrameLabel->{label[[paper1]],"Probability",Distr[[paper1]],
                  "Return Period"},PlotJoined->True,
              PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
              DisplayFunction:>Identity];
          
          DD9=ListPlot[Bandminus,
              FrameTicks->{Grid1,Probabilities1,Automatic,ReturnPeriods},
              Frame->True,
              GridLines->{Grid2,PTran[paper1,Probabilities]},
              Prolog->AbsolutePointSize[4],
              FrameLabel->{label[[paper1]],"Probability",Distr[[paper1]],
                  "Return Period"},PlotJoined->True,
              PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
              DisplayFunction:>Identity];
          
          DD1=ListPlot[ExtremesLine,
              FrameTicks->{Grid1,Probabilities1,Automatic,ReturnPeriods},
              Frame->True,
              GridLines->{Grid2,PTran[paper1,Probabilities]},
              Prolog->AbsolutePointSize[4],
              FrameLabel->{label[[paper1]],"Probability",Distr[[paper1]],
                  "Return Period"},PlotJoined->True,
              PlotStyle->AbsoluteThickness[1],
              DisplayFunction:>Identity];
          
          DD2=ListPlot[Points,
              FrameTicks->{Grid1,Probabilities1,Automatic,ReturnPeriods},
              GridLines->{Grid2,PTran[paper1,Probabilities]},
              Prolog->AbsolutePointSize[4],
              FrameLabel->{label[[paper1]],"Probability",Distr[[paper1]],
                  "Return Period"},DisplayFunction:>Identity],
          DD8=ListPlot[Bandplus,
              FrameTicks->{Grid1,Probabilities1,Automatic,ReturnPeriods},
              AxesOrigin->{Points[[1,1]],Points[[1,2]]},Frame->True,
              GridLines->{Grid2,PTran[paper1,Probabilities]},
              Prolog->AbsolutePointSize[4],
              FrameLabel->{"X","Probability",Distr[[paper1]],
                  "Return Period"},PlotJoined->True,
              PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
              DisplayFunction:>Identity];
          
          DD9=ListPlot[Bandminus,
              FrameTicks->{Grid1,Probabilities1,Automatic,ReturnPeriods},
              AxesOrigin->{Points[[1,1]],Points[[1,2]]},Frame->True,
              GridLines->{Grid2,PTran[paper1,Probabilities]},
              Prolog->AbsolutePointSize[4],
              FrameLabel->{"X","Probability",Distr[[paper1]],
                  "Return Period"},PlotJoined->True,
              PlotStyle->{AbsoluteThickness[1],Dashing[{0.02,0.02}]},
              DisplayFunction:>Identity];
          
          DD1=ListPlot[ExtremesLine,
              FrameTicks->{Grid1,Probabilities1,Automatic,ReturnPeriods},
              AxesOrigin->{Points[[1,1]],Points[[1,2]]},Frame->True,
              GridLines->{Grid2,PTran[paper1,Probabilities]},
              Prolog->AbsolutePointSize[4],
              FrameLabel->{"X","Probability",Distr[[paper1]],
                  "Return Period"},PlotJoined->True,
              PlotStyle->AbsoluteThickness[1],
              DisplayFunction:>Identity];
          
          DD2=ListPlot[Points,
              PlotRange->{{Points[[1,1]],
                    Points[[1,2]]},{Points[[SampleSize,1]],
                    Points[[SampleSize,2]]}},
              FrameTicks->{Grid1,Probabilities1,Automatic,ReturnPeriods},
              AxesOrigin->{Points[[1,1]],Points[[1,2]]},
              GridLines->{Grid2,PTran[paper1,Probabilities]},
              Prolog->AbsolutePointSize[4],
              FrameLabel->{"X","Probability",Distr[[paper1]],
                  "Return Period"},DisplayFunction:>Identity]],
        Return[]];
      DD3=
        Show[DD1,DD2,DD8,DD9,Axes->False,Frame->True,
          DisplayFunction:>$DisplayFunction,PlotRange->All,
          ImageSize->{imSize,Automatic}];
      Return[DD3]];
         
         
ASAEGEV[Sample_,lambda_,delta_,kappa_,Tail_]:=
  Module[ {aux,i,Sample1,n=Length[Sample]},
  If[Tail==1,Sample1=Sort[Sample],Sample1=-Sort[Sample]];
    If[kappa!=0,
      aux=Sum[Abs[(Sample1[[i]]-lambda-
                    delta*(1-(-Log[(i-0.5)/n])^kappa)/kappa)]/(Sample1[[n]]-
                  Sample1[[1]]),{i,1,n}]/n,
      Sum[Abs[(Sample1[[i]]-lambda+delta*Log[-Log[(i-0.5)/n]])]/(Sample1[[n]]-\

                Sample1[[1]]),{i,1,n}]/n];Return[aux]
    ];
    
      ASAEGPD[Sample_,lambda_,kappa_,Tail_]:=
  Module[{Sample1,i,aux=0,n=Length[Sample]},
    If[Tail==1,Sample1=Sort[Sample],Sample1=-Sort[Sample]];
    If[Abs[kappa]>0.01,
      aux=Sum[Abs[Sample1[[i]]-lambda*(1-(1-(i-0.5)/n)^kappa)/kappa],{i,1,
              n}]/(n*(Sample1[[n]]-Sample1[[1]])),
      aux=Sum[Abs[Sample1[[i]]+lambda*Log[1-(i-0.5)/n]],{i,1,
              n}]/(n*(Sample1[[n]]-Sample1[[1]]))];
    Return[aux]
    ];


    
(*****************************************************************************\


*****)
 

    
End[]

EndPackage[]

