{"word_count": 728, "character_count": 7256}
These data were collected at Copenhagen Reinsurance and comprise 2167 fire losses over the period 1980 to 1990, They have been adjusted for inflation to reflect 1985 values and are expressed in millions of Danish Kron. Note that it is possible to work with the same data as above but the total claim has been divided into a building loss, a loss of contents and a loss of profits.&#62; base1=read.table(+ "http://freakonometrics.free.fr/danish-univariate.txt",+ header=TRUE)&#62; base2=read.table(+ "http://freakonometrics.free.fr/danish-multivariate.txt",+ header=TRUE)Consider here the first dataset (we deal - so far - with univariate extremes),&#62; X=base1$Loss.in.DKM&#62; D=as.Date(as.character(base1$Date),"%m/%d/%Y")&#62; plot(D,X,type="h")The graph is the following,A natural idea is then to ploti.e.&#62; Xs=sort(X)&#62; logXs=rev(log(Xs))&#62; n=length(X)&#62; plot(log(Xs),log((n:1)/(n+1)))Points are on a straight line here. The slope can be obtained using a linear regression,&#62; B=data.frame(X=log(Xs),Y=log((n:1)/(n+1)))&#62; reg=lm(Y~X,data=B)&#62; summary(reg)Call:lm(formula = Y ~ X, data = B)Residuals:Min       1Q   Median       3Q      Max-0.59999 -0.00777  0.00878  0.02461  0.20309Coefficients:Estimate Std. Error t value Pr(&#62;|t|)(Intercept)  0.089442   0.001572   56.88   &#60;2e-16 ***X           -1.382181   0.001477 -935.55   &#60;2e-16 ***---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.04928 on 2165 degrees of freedomMultiple R-squared: 0.9975,Adjusted R-squared: 0.9975F-statistic: 8.753e+05 on 1 and 2165 DF,  p-value: &#60; 2.2e-16&#62; reg=lm(Y~X,data=B[(n-500):n,])&#62; summary(reg)Call:lm(formula = Y ~ X, data = B[(n - 500):n, ])Residuals:Min       1Q   Median       3Q      Max-0.48502 -0.02148 -0.00900  0.01626  0.35798Coefficients:Estimate Std. Error t value Pr(&#62;|t|)(Intercept)  0.186188   0.010033   18.56   &#60;2e-16 ***X           -1.432767   0.005105 -280.68   &#60;2e-16 ***---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.07751 on 499 degrees of freedomMultiple R-squared: 0.9937,Adjusted R-squared: 0.9937F-statistic: 7.878e+04 on 1 and 499 DF,  p-value: &#60; 2.2e-16&#62; reg=lm(Y~X,data=B[(n-100):n,])&#62; summary(reg)Call:lm(formula = Y ~ X, data = B[(n - 100):n, ])Residuals:Min       1Q   Median       3Q      Max-0.33396 -0.03743  0.02279  0.04754  0.62946Coefficients:Estimate Std. Error t value Pr(&#62;|t|)(Intercept)  0.67377    0.06777   9.942   &#60;2e-16 ***X           -1.58536    0.02240 -70.772   &#60;2e-16 ***---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.1299 on 99 degrees of freedomMultiple R-squared: 0.9806,Adjusted R-squared: 0.9804F-statistic:  5009 on 1 and 99 DF,  p-value: &#60; 2.2e-16The slope here is somehow related to the tail index of the distribution. Consider some heavy tailed distribution, i.e. , so that , where  is some slowly varying function. Equivalently, the exists a slowly varying function  such that . Theni.e. since a natural estimator for  is the order statistic , the slope of the straight line is the opposite of tail index . The estimator of the slope is (considering only the  largest observations)Hill's estimator is based on the assumption that the denominator above is almost 1 (which means that  , as ), i.e.Note that, if , but not two fast, i.e.  as , then  (one can even get   with stronger convergence assumptions). FurtherBased on that (asymptotic) distribution, it is possible to get a (asymptotic) confidence interval for &#62; xi=1/(1:n)*cumsum(logXs)-logXs&#62; xise=1.96/sqrt(1:n)*xi&#62; plot(1:n,xi,type="l",ylim=range(c(xi+xise,xi-xise)),+ xlab="",ylab="",)&#62; polygon(c(1:n,n:1),c(xi+xise,rev(xi-xise)),+ border=NA,col="lightblue")&#62; lines(1:n,xi+xise,col="red",lwd=1.5)&#62; lines(1:n,xi-xise,col="red",lwd=1.5)&#62; lines(1:n,xi,lwd=1.5)&#62; abline(h=0,col="grey")It is also possible to work with , then . And similarly  as  (and again  with additional assumptions on the rate of convergence), and(obtained using the delta-method). Again, we can use that result to derive (asymptotic) confidence intervals&#62; alpha=1/xi&#62; alphase=1.96/sqrt(1:n)/xi&#62; YL=c(0,3)&#62; plot(1:n,alpha,type="l",ylim=YL,xlab="",ylab="",)&#62; polygon(c(1:n,n:1),c(alpha+alphase,rev(alpha-alphase)),+ border=NA,col="lightblue")&#62; lines(1:n,alpha+alphase,col="red",lwd=1.5)&#62; lines(1:n,alpha-alphase,col="red",lwd=1.5)&#62; lines(1:n,alpha,lwd=1.5)&#62; abline(h=0,col="grey")The Deckers-Einmahl-de Haan estimator iswhere forThen (given again conditions on the speed of convergence i.e. , with  as ),Finally, Pickands' estimatorit is possible to prove that, as ,Here the code is&#62; Xs=rev(sort(X))&#62; xi=1/log(2)*log( (Xs[seq(1,length=trunc(n/4),by=1)]-+ Xs[seq(2,length=trunc(n/4),by=2)])/+ (Xs[seq(2,length=trunc(n/4),by=2)]-Xs[seq(4,+ length=trunc(n/4),by=4)]) )&#62; xise=1.96/sqrt(seq(1,length=trunc(n/4),by=1))*+sqrt( xi^2*(2^(xi+1)+1)/((2*(2^xi-1)*log(2))^2))&#62; plot(seq(1,length=trunc(n/4),by=1),xi,type="l",+ ylim=c(0,3),xlab="",ylab="",)&#62; polygon(c(seq(1,length=trunc(n/4),by=1),rev(seq(1,+ length=trunc(n/4),by=1))),c(xi+xise,rev(xi-xise)),+ border=NA,col="lightblue")&#62; lines(seq(1,length=trunc(n/4),by=1),+ xi+xise,col="red",lwd=1.5)&#62; lines(seq(1,length=trunc(n/4),by=1),+ xi-xise,col="red",lwd=1.5)&#62; lines(seq(1,length=trunc(n/4),by=1),xi,lwd=1.5)&#62; abline(h=0,col="grey")It is also possible to use maximum likelihood techniques to fit a GPD distribution over a high threshold.&#62; library(evd)&#62; library(evir)&#62; gpd(X,5)$n[1] 2167$threshold[1] 5$p.less.thresh[1] 0.8827873$n.exceed[1] 254$method[1] "ml"$par.estsxi      beta0.6320499 3.8074817$par.sesxi      beta0.1117143 0.4637270$varcov[,1]        [,2][1,]  0.01248007 -0.03203283[2,] -0.03203283  0.21504269$information[1] "observed"$converged[1] 0$nllh.final[1] 754.1115attr(,"class")[1] "gpd"or equivalently (or almost)&#62; gpd.fit(X,5)$threshold[1] 5$nexc[1] 254$conv[1] 0$nllh[1] 754.1115$mle[1] 3.8078632 0.6315749$rate[1] 0.1172127$se[1] 0.4636270 0.1116136The interest of the latest function is that it is possible to visualize the profile likelihood of the tail index,&#62; gpd.profxi(gpd.fit(X,5),xlow=0,xup=3)or&#62; gpd.profxi(gpd.fit(X,20),xlow=0,xup=3)Hence, it is possible to plot the maximum likelihood estimator of the tail index, as a function of the threshold (including a confidence interval),&#62; GPDE=Vectorize(function(u){gpd(X,u)$par.ests[1]})&#62; GPDS=Vectorize(function(u){+ gpd(X,u)$par.ses[1]})&#62; u=c(seq(2,10,by=.5),seq(11,25))&#62; XI=GPDE(u)&#62; XIS=GPDS(u)&#62; plot(u,XI,ylim=c(0,2))&#62; segments(u,XI-1.96*XIS,u,XI++ 1.96*XIS,lwd=2,col="red")Finally, it is possible to use block-maxima techniques.&#62; gev.fit(X)$conv[1] 0$nllh[1] 3392.418$mle[1] 1.4833484 0.5930190 0.9168128$se[1] 0.01507776 0.01866719 0.03035380The estimator of the tail index is here the last coefficient, on the right.Since it is rather difficult to install a package in class rooms, here is the source of rcodes used here (to fit a GPD for exceedances)&#62; source("http://freakonometrics.blog.free.fr/public/code/gpd.R")Next time, we will discuss how to use those estimators.