1 Parties intermédiaires

Début On peut maintenant retracer le qqplot avec les seuils choisis et examiner si les valeurs entre les deux seuils sont normalement distribuées :

probs=seq(0.0001,0.9999,1e-4)

qqd=data.frame(probs=probs,
               th=qnorm(p=probs,
                        mean=mean(data$DR[data$DR>-2.19 & data$DR<2.06],na.rm=TRUE),
                        sd=sd(data$DR[data$DR>-2.19 & data$DR<2.06],na.rm=TRUE)),
               emp=quantile(data$DR,probs,na.rm=TRUE))

p <- ggplot(qqd, aes(th,emp,label=probs))
p=p+geom_point()+geom_abline(slope=1,intercept = 0,colour=2)

ggplotly(p)

On voit qu’il est nécessaire de choisir un autre type de lois pour les parties intermédiaires qui sont entre les valeurs extrêmes et les valeurs normalement distribuées.

2 Modélisation des queues de distributions

2.1 Loi des valeurs extrêmes

2.1.1 Loi des rendements postifs extrêmes

seuilp2=2.06
rp=data$DR[data$DR>seuilp2 & is.na(data$DR)==FALSE]
length(rp)
## [1] 323
fitg=fitdist(rp-seuilp2,method="mle",distr="gamma")
fitw=fitdist(rp-seuilp2,method="mle",distr="weibull")

probs=seq(0.0001,0.9999,1/333)

qqd=rbind(data.frame(probs=probs,
                     emp=quantile(rp,probs,na.rm=TRUE),
                     th=qpareto(p=probs,
                        scale=2.06,
                        shape=3.3),
               loi="pareto"),
          data.frame(probs=probs,
                     emp=quantile(rp,probs,na.rm=TRUE),
                     th=qgamma(p=probs,
                        shape=fitg$estimate[1],
                        rate=fitg$estimate[2])+seuilp2,
               loi="gamma"),
          data.frame(probs=probs,
                     emp=quantile(rp,probs,na.rm=TRUE),
                     th=qweibull(p=probs,
                        shape=fitw$estimate[1],
                        scale=fitw$estimate[2])+seuilp2,
               loi="weibull"))

p <- ggplot(data=qqd, aes(th,emp,label=probs,colour=loi))
p=p+geom_point()+
  geom_abline(slope=1,intercept = 0,colour=2)

ggplotly(p)

2.1.2 Loi des rendements négatifs extrêmes

Même chose, non affiché pour alléger le fichier.

2.2 Modélisation des fréquences

On étudie la fréquence annuelle des rendements extrêmes.

freq=function(data,seuil1,seuil2){
  d=data[data$DR>seuil1 & data$DR<seuil2 & is.na(data$DR)==FALSE,]
  freq=merge(d[,list(nb=.N),by="Year"],data[,list(nb=.N),by="Year"],by="Year")
  freq=freq$nb.x/freq$nb.y
  mf=mean(freq)*252
  vf=var(freq)*252^2
  sdf=sqrt(vf)
  return( list(mf=mf,vf=vf,sdf=sdf))
}

fep=freq(data,2.06,Inf)
fen=freq(data,-Inf,-2.06)
fep
## $mf
## [1] 6.635477
## 
## $vf
## [1] 46.37949
## 
## $sdf
## [1] 6.810249
## Simulation des nombres par intervalle
OEP=function(seuil,alpha,mf,vf){
  layers=seq(30,45,0.5)
  ns=100000
  res=matrix(0,ns,2)
  nlayers=length(layers)-1
  nb_m = matrix(0,ns,nlayers)
  
  for (k in 1:ns){
    nb=rnbinom(1,size=mf^2/(vf-mf),mu=mf)
    te=rpareto(nb,seuil,alpha)
    mx=max(te)
    res[k,1]=nb
    res[k,2]=mx
    nb_t=NULL
    for (i in 1:(nlayers)) {
      nb_t = c(nb_t,sum((te < layers[i + 1])*(te >= layers[i])))
    }
    nb_m[k,]=nb_t
  }
  nb_m[is.na(nb_m)==TRUE]=0
  layer_freq = colMeans(nb_m)
  layer_freq
  rev(cumsum(rev(layer_freq)))
  return (list(res=res,layer_freq=layer_freq))
}

Copyright © 2016 Blog de Kezhan Shi