This Tuesday, I will be giving the second part of the (crash) graduate course on advanced tools for econometrics. It will take place in Rennes, IMAPP room, and I have been told that there will be a visio with Nantes and Angers. Slides for the morning are online, as well as slides for the afternoon.
In the morning, we will talk about variable section and penalization, and in the afternoon, it will be on changing the loss function (quantile regression).
Exponential Smoothing is an old technique, but it can perform extremely well on real time series, as discussed in Hyndman, Koehler, Ord & Snyder (2008)),
when Gardner (2005) appeared, many believed that exponential smoothing should be disregarded because it was either a special case of ARIMA modeling or an ad hoc procedure with no statistical rationale. As McKenzie (1985) observed, this opinion was expressed in numerous references to my paper. Since 1985, the special case argument has been turned on its head, and today we know that exponential smoothing methods are optimal for a very general class of state-space models that is in fact broader than the ARIMA class.
Furthermore, I like it because I think it has nice pedagogical features. Consider simple exponential smoothing, L_{t}=\alpha Y_{t}+(1-\alpha)L_{t-1} where \alpha\in(0,1) is the smoothing weight. It is locally constant, in the sense that {}_{t}\hat Y_{t+h} = L_{t}
library(datasets) X=as.numeric(Nile) SimpleSmooth = function(a){ T=length(X) L=rep(NA,T) L[1]=X[1] for(t in 2:T){L[t]=a*X[t]+(1-a)*L[t-1]} return(L) } plot(X,type="b",cex=.6) lines(SimpleSmooth(.2),col="red") |
When using the standard R function, we get
hw=HoltWinters(X,beta=FALSE,gamma=FALSE, l.start=X[1]) hw$alpha [1] 0.2465579 |
Of course, one can replicate that optimal value
V=function(a){ T=length(X) L=erreur=rep(NA,T) erreur[1]=0 L[1]=X[1] for(t in 2:T){ L[t]=a*X[t]+(1-a)*L[t-1] erreur[t]=X[t]-L[t-1] } return(sum(erreur^2)) } optim(.5,V)$par [1] 0.2464844 |
Here, the optimal value for \alpha is the one that minimizes the one-step prediction, for the \ell_2 loss function, i.e. \sum_{t=2}^n(Y_t-{}_{t-1}\hat Y_t)^2 where here {}_{t-1}\hat Y_t = L_{t-1}. But one can consider another loss function, for instance the quantile loss function, \ell_{\tau}(\varepsilon)=\varepsilon(\tau-\mathbb{I}_{\varepsilon\leq 0}). The optimal coefficient is then obtained using
HWtau=function(tau){ loss=function(e) e*(tau-(e<=0)*1) V=function(a){ T=length(X) L=erreur=rep(NA,T) erreur[1]=0 L[1]=X[1] for(t in 2:T){ L[t]=a*X[t]+(1-a)*L[t-1] erreur[t]=X[t]-L[t-1] } return(sum(loss(erreur))) } optim(.5,V)$par } |
Here is the evolution of \alpha^\star_\tau as a function of \tau (the level of the quantile considered).
T=(1:49)/50 HW=Vectorize(HWtau)(T) plot(T,HW,type="l") abline(h= hw$alpha,lty=2,col="red") |
Note that the optimal \alpha is decreasing with \tau. I wonder how general this result can be…
Of course, one can consider more general exponential smoothing, for instance the double one, with L_t=\alpha Y_t+(1-\alpha)[L_{t-1}+B_{t-1}]andB_t=\beta[L_t-L_{t-1}]+(1-\beta)B_{t-1}so that the prediction is now {}_{t}\hat Y_{t+h} = L_{t}+hB_t (it is now locally linear – and no longer constant).
hw=HoltWinters(X,gamma=FALSE,l.start=X[1]) hw$alpha alpha 0.4200241 hw$beta beta 0.05973389 |
The code to compute the smoothed series is the following
DoubleSmooth = function(a,b){ T=length(X) L=B=rep(NA,T) L[1]=X[1]; B[1]=0 for(t in 2:T){ L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1]) B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] } return(L+B) } |
Here also it is possible to replicate R using the \ell_2 loss function
V=function(A){ a=A[1] b=A[2] T=length(X) L=B=erreur=rep(NA,T) erreur[1]=0 L[1]=X[1]; B[1]=X[2]-X[1] for(t in 2:T){ L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1]) B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] erreur[t]=X[t]-(L[t-1]+B[t-1]) } return(sum(erreur^2)) } optim(c(.5,.05),V)$par [1] 0.41904510 0.05988304 |
(up to numerical optimization approximation, I guess). But here also, a quantile loss function can be considered
HWtau=function(tau){ loss=function(e) e*(tau-(e<=0)*1) V=function(A){ a=A[1] b=A[2] T=length(X) L=B=erreur=rep(NA,T) erreur[1]=0 L[1]=X[1]; B[1]=X[2]-X[1] for(t in 2:T){ L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1]) B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] erreur[t]=X[t]-(L[t-1]+B[t-1]) } return(sum(loss(erreur))) } optim(c(.5,.05),V)$par } |
and we can plot those values on a graph
T=(1:49)/50 HW=Vectorize(HWtau)(T) plot(HW[1,],HW[2,],type="l") abline(v= hw$alpha,lwd=.4,lty=2,col="red") abline(h= hw$beta,lwd=.4,lty=2,col="red") points(hw$alpha,hw$beta,pch=19,col="red") |
(with \alpha on the x-axis, and \beta on the y-axis). So here, it is extremely simple to change the loss function, but so far, it should be done manually. Of course, one do it also for the seasonal exponential smoothing model.
Consider our loss-ALAE dataset, and – as in Frees & Valdez (1998) – let us fit a parametric model, in order to price a reinsurance treaty. The dataset is the following,
> library(evd) > data(lossalae) > Z=lossalae > X=Z[,1];Y=Z[,2]
The first step can be to estimate marginal distributions, independently. Here, we consider lognormal distributions for both components,
> Fempx=function(x) mean(X<=x) > Fx=Vectorize(Fempx) > u=exp(seq(2,15,by=.05)) > plot(u,Fx(u),log="x",type="l", + xlab="loss (log scale)") > Lx=function(px) -sum(log(Vectorize(dlnorm)( + X,px[1],px[2]))) > opx=optim(c(1,5),fn=Lx) > opx$par [1] 9.373679 1.637499 > lines(u,Vectorize(plnorm)(u,opx$par[1], + opx$par[2]),col="red")
The fit here is quite good,
For the second component, we do the same,
> Fempy=function(x) mean(Y<=x) > Fy=Vectorize(Fempy) > u=exp(seq(2,15,by=.05)) > plot(u,Fy(u),log="x",type="l", + xlab="ALAE (log scale)") > Ly=function(px) -sum(log(Vectorize(dlnorm)( + Y,px[1],px[2]))) > opy=optim(c(1.5,10),fn=Ly) > opy$par [1] 8.522452 1.429645 > lines(u,Vectorize(plnorm)(u,opy$par[1], + opy$par[2]),col="blue")
It is not as good as the fit obtained on losses, but it is not that bad,
Now, consider a multivariate model, with Gumbel copula. We’ve seen before that it worked well. But this time, consider the maximum likelihood estimator globally.
> Cop=function(u,v,a) exp(-((-log(u))^a+ + (-log(v))^a)^(1/a)) > phi=function(t,a) (-log(t))^a > cop=function(u,v,a) Cop(u,v,a)*(phi(u,a)+ + phi(v,a))^(1/a-2)*( + a-1+(phi(u,a)+phi(v,a))^(1/a))*(phi(u,a-1)* + phi(v,a-1))/(u*v) > L=function(p) {-sum(log(Vectorize(dlnorm)( + X,p[1],p[2])))- + sum(log(Vectorize(dlnorm)(Y,p[3],p[4])))- + sum(log(Vectorize(cop)(plnorm(X,p[1],p[2]), + plnorm(Y,p[3],p[4]),p[5])))} > opz=optim(c(1.5,10,1.5,10,1.5),fn=L) > opz$par [1] 9.377219 1.671410 8.524221 1.428552 1.468238
Marginal parameters are (slightly) different from the one obtained independently,
> c(opx$par,opy$par) [1] 9.373679 1.637499 8.522452 1.429645 > opz$par[1:4] [1] 9.377219 1.671410 8.524221 1.428552
And the parameter of Gumbel copula is close to the one obtained with heuristic methods in class.
Now that we have a model, let us play with it, to price a reinsurance treaty. But first, let us see how to generate Gumbel copula… One idea can be to use the frailty approach, based on a stable frailty. And we can use Chambers et al (1976)to generate a stable distribution. So here is the algorithm to generate samples from Gumbel copula
> alpha=opz$par[5] > invphi=function(t,a) exp(-t^(1/a)) > n=500 > x=matrix(rexp(2*n),n,2) > angle=runif(n,0,pi) > E=rexp(n) > beta=1/alpha > stable=sin((1-beta)*angle)^((1-beta)/beta)* + (sin(beta*angle))/(sin(angle))^(1/beta)/ + (E^(alpha-1)) > U=invphi(x/stable,alpha) > plot(U)
Here, we consider only 500 simulations,
> Xloss=qlnorm(U[,1],opz$par[1],opz$par[2]) > Xalae=qlnorm(U[,2],opz$par[3],opz$par[4])
In standard reinsurance treaties – see e.g. Clarke (1996) – allocated expenses are splited prorata capita between the insurance company, and the reinsurer. If denotes losses, and
the allocated expenses, a standard excess treaty can be has payoff
> L=100000 > R=50000 > Z=((Xloss-R)+(Xloss-R)/Xloss*Xalae)* + (R<=Xloss)*(Xloss<L)+ + ((L-R)+(L-R)/R*Xalae)*(L<=Xloss) > mean(Z) [1] 12596.45
Now, play with it… it is possible to find a better fit, I guess…
Today, we will go further on the inference of copula functions. Some codes (and references) can be found on a previous post, on nonparametric estimators of copula densities (among other related things). Consider (as before) the loss-ALAE dataset (since we’ve been working a lot on that dataset)
> library(MASS) > library(evd) > X=lossalae > U=cbind(rank(X[,1])/(nrow(X)+1),rank(X[,2])/(nrow(X)+1))
The standard tool to plot nonparametric estimators of densities is to use multivariate kernels. We can look at the density using
> mat1=kde2d(U[,1],U[,2],n=35) > persp(mat1$x,mat1$y,mat1$z,col="green", + shade=TRUE,theta=s*5, + xlab="",ylab="",zlab="",zlim=c(0,7))
or level curves (isodensity curves) with more detailed estimators (on grids with shorter steps)
> mat1=kde2d(U[,1],U[,2],n=101) > image(mat1$x,mat1$y,mat1$z,col= + rev(heat.colors(100)),xlab="",ylab="") > contour(mat1$x,mat1$y,mat1$z,add= + TRUE,levels = pretty(c(0,4), 11))
Kernels are nice, but we clearly observe some border bias, extremely strong in corners (the estimator is 1/4th of what it should be, see another post for more details). Instead of working on sample on the unit square, consider some transformed sample
, where
is a given function. E.g. a quantile function of an unbounded distribution, for instance the quantile function of the
distribution. Then, we can estimate the density of the transformed sample, and using the inversion technique, derive an estimator of the density of the initial sample. Since the inverse of a (general) function is not that simple to compute, the code might be a bit slow. But it does work,
> gaussian.kernel.copula.surface <- function (u,v,n) { + s=seq(1/(n+1), length=n, by=1/(n+1)) + mat=matrix(NA,nrow = n, ncol = n) + sur=kde2d(qnorm(u),qnorm(v),n=1000, + lims = c(-4, 4, -4, 4)) + su<-sur$z + for (i in 1:n) { + for (j in 1:n) { + Xi<-round((qnorm(s[i])+4)*1000/8)+1; + Yj<-round((qnorm(s[j])+4)*1000/8)+1 + mat[i,j]<-su[Xi,Yj]/(dnorm(qnorm(s[i]))* + dnorm(qnorm(s[j]))) + } + } + return(list(x=s,y=s,z=data.matrix(mat))) + }
Here, we get
Note that it is possible to consider another transformation, e.g. the quantile function of a Student-t distribution.
> student.kernel.copula.surface = + function (u,v,n,d=4) { + s <- seq(1/(n+1), length=n, by=1/(n+1)) + mat <- matrix(NA,nrow = n, ncol = n) + sur<-kde2d(qt(u,df=d),qt(v,df=d),n=5000, + lims = c(-8, 8, -8, 8)) + su<-sur$z + for (i in 1:n) { + for (j in 1:n) { + Xi<-round((qt(s[i],df=d)+8)*5000/16)+1; + Yj<-round((qt(s[j],df=d)+8)*5000/16)+1 + mat[i,j]<-su[Xi,Yj]/(dt(qt(s[i],df=d),df=d)* + dt(qt(s[j],df=d),df=d)) + } + } + return(list(x=s,y=s,z=data.matrix(mat))) + }
Another strategy is to consider kernel that have precisely the unit interval as support. The idea is here to consider the product of Beta kernels, where parameters depend on the location
> beta.kernel.copula.surface= + function (u,v,bx=.025,by=.025,n) { + s <- seq(1/(n+1), length=n, by=1/(n+1)) + mat <- matrix(0,nrow = n, ncol = n) + for (i in 1:n) { + a <- s[i] + for (j in 1:n) { + b <- s[j] + mat[i,j] <- sum(dbeta(a,u/bx,(1-u)/bx) * + dbeta(b,v/by,(1-v)/by)) / length(u) + } + } + return(list(x=s,y=s,z=data.matrix(mat))) + }
On those two graphs, we can clearly observe strong tail dependence in the upper (right) corner, that cannot be intuited using a standard kernel estimator…
As mentioned in the course last week Venter (2003) suggested nice functions to illustrate tail dependence (see also some slides used in Berlin a few years ago).
Joe (1990) suggested a (strong) tail dependence index. For lower tails, for instance, consider
i.e
The idea is to plot the function above, in order to visualize limiting behavior. Define
for the lower tail, and
is the survival copula associated with 
, in the sense that
while
Now, one can easily derive empirical conterparts of those function, i.e.
and
Thus, for upper tail, on the right, we have the following graph
and for the lower tail, on the left, we have
For the code, consider some real data, like the loss-ALAE dataset.
> library(evd) > X=lossalae
The idea is to plot, on the left, the lower tail concentration function, and on the right, the upper tail function.
> U=rank(X[,1])/(nrow(X)+1) > V=rank(X[,2])/(nrow(X)+1) > Lemp=function(z) sum((U<=z)&(V<=z))/sum(U<=z) > Remp=function(z) sum((U>=1-z)&(V>=1-z))/sum(U>=1-z) > u=seq(.001,.5,by=.001) > L=Vectorize(Lemp)(u) > R=Vectorize(Remp)(rev(u)) > plot(c(u,u+.5-u[1]),c(L,R),type="l",ylim=0:1, + xlab="LOWER TAIL UPPER TAIL") > abline(v=.5,col="grey")
Now, we can compare this graph, with what should be obtained for some parametric copulas that have the same Kendall’s tau (e.g.). For instance, if we consider a Gaussian copula,
> tau=cor(lossalae,method="kendall")[1,2] > library(copula) > paramgauss=sin(tau*pi/2) > copgauss=normalCopula(paramgauss) > Lgaussian=function(z) pCopula(c(z,z),copgauss)/z > Rgaussian=function(z) (1-2*z+pCopula(c(z,z),copgauss))/(1-z) > u=seq(.001,.5,by=.001) > Lgs=Vectorize(Lgaussian)(u) > Rgs=Vectorize(Rgaussian)(1-rev(u)) > lines(c(u,u+.5-u[1]),c(Lgs,Rgs),col="red")
or Gumbel’s copula,
> paramgumbel=1/(1-tau) > copgumbel=gumbelCopula(paramgumbel, dim = 2) > Lgumbel=function(z) pCopula(c(z,z),copgumbel)/z > Rgumbel=function(z) (1-2*z+pCopula(c(z,z),copgumbel))/(1-z) > u=seq(.001,.5,by=.001) > Lgl=Vectorize(Lgumbel)(u) > Rgl=Vectorize(Rgumbel)(1-rev(u)) > lines(c(u,u+.5-u[1]),c(Lgl,Rgl),col="blue")
That’s nice (isn’t it?), but since we do not have any confidence interval, it is still hard to conclude (even if it looks like Gumbel copula has a much better fit than the Gaussian one). A strategy can be to generate samples from those copulas, and to visualize what we had. With a Gaussian copula, the graph looks like
> u=seq(.0025,.5,by=.0025); nu=length(u) > nsimul=500 > MGS=matrix(NA,nsimul,2*nu) > for(s in 1:nsimul){ + Xs=rCopula(nrow(X),copgauss) + Us=rank(Xs[,1])/(nrow(Xs)+1) + Vs=rank(Xs[,2])/(nrow(Xs)+1) + Lemp=function(z) sum((Us<=z)&(Vs<=z))/sum(Us<=z) + Remp=function(z) sum((Us>=1-z)&(Vs>=1-z))/sum(Us>=1-z) + MGS[s,1:nu]=Vectorize(Lemp)(u) + MGS[s,(nu+1):(2*nu)]=Vectorize(Remp)(rev(u)) + lines(c(u,u+.5-u[1]),MGS[s,],col="red") + }
(including – pointwise – 90% confidence bands)
> Q95=function(x) quantile(x,.95) > V95=apply(MGS,2,Q95) > lines(c(u,u+.5-u[1]),V95,col="red",lwd=2) > Q05=function(x) quantile(x,.05) > V05=apply(MGS,2,Q05) > lines(c(u,u+.5-u[1]),V05,col="red",lwd=2)
while it is
with Gumbel copula. Isn’t it a nice (graphical) tool ?
But as mentioned in the course, the statistical convergence can be slow. Extremely slow. So assessing if the underlying copula has tail dependence, or not, it now that simple. Especially if the copula exhibits tail independence. Like the Gaussian copula. Consider a sample of size 1,000. This is what we obtain if we generate random scenarios,
or we look at the left tail (with a log-scale)
Now, consider a 10,000 sample,
or with a log-scale
We can even consider a 100,000 sample,
or with a log-scale
On those graphs, it is rather difficult to conclude if the limit is 0, or some strictly positive value (again, it is a classical statistical problem when the value of interest is at the border of the support of the parameter). So, a natural idea is to consider a weaker tail dependence index. Unless you have something like 100,000 observations…
Suite à une formation que je faisais en fin de semaine à Brest (les slides sont ici et là), je voulais revenir sur les histoires de tails of copulas, pour reprendre le titre de l’article (ici) de Gary Venter (et qui correspond à des choses que j’avais pu présenter il y a quelques années à Berlin, les slides étant en ligne ici).
L’idée est de noter qu’il est noter qu’il existe deux manières de quantifier la dépendance de queue. La première est liée à l’approche de Joe (1990, ici, ou 1997 pour le livre), qui a introduit un (strong) tail dependence index. Par exemple pour la queue inférieure,
soit
La seconde est liée à une idée que l’on retrouve dans les travaux de Janet Heffernan, Stuart Coles ou Jonathan Tawn. L’intuition est la suivante (on peut la retrouver en ligne ici). Si 
et 
ont la même loi et que l’on suppose les variables indépendantes, alors
En revanche, si les variables sont comonotones (c’est à dire égales comme on suppose les lois identiques),
Aussi, on peut supposer qu’il existe un indice 
tel que
Le soucis est que le cas d’indépendance correspond à =2, alors que le cas de dépendance forte correspond au cas =1. Il est alors usuel de faire une transformation affine pour se ramener sur [0,1], et que la force de la dépendance soit croissante avec , e.g.
Posons alors
qui pourra être interprété comme un (weak) tail dependence index.
Bref, ces deux mesures donnent de l’information sur le comportement dans les queues de distribution.
L’idée est de noter qu’il est possible d’étudier ces fonctions afin de mieux comprendre le comportement dans les queues. En s’inspirant de Gary Venter, on peut définir
pour étudier le comportement dans la queue inférieure, et
est la copule de survie associée à 
, au sens où
et
Cet outil permettra de modéliser la dépendance forte. On peut également poser, afin d’étudier la dépendance faible,
ou
L’idée est de noter qu’il est facile d’estimer ces fonctions. Ces outils peuvent être utiles pour mieux comprendre le comportement dans les queues.
Par exemple pour une copule Gaussienne de corrélation 0,5, on a la forme théorique suivante pour les fonctions de concentration (au sens fort)
Statistiquement, il est possible d’estimer ces quantités en comptant simplement le nombre d’observations dans le coin inférieur gauche, ou le coin supérieur droit. Si on dispose d’un échantillon, on peut alors regarder ce que donnent les versions
et
Pour un échantillon de taille n=500, on obtient les intervalles de confiance à 90% de la forme suivante,
Le code R ressemble à ça
> library(evd); data(lossalae)
> cor(lossalae,method="spearman")
Loss ALAE
Loss 1.000000 0.451872
ALAE 0.451872 1.000000
avec le code suivant pour la version empirique,
> z=seq(0,.5,by=.001) > U=rank(v[,1])/(nrow(v)+1) > V=rank(v[,2])/(nrow(v)+1) > Lemp=rep(NA,length(z)) > Remp=rep(NA,length(z)) > for(i in 1:length(z)){ + Lemp[i]=sum((U<=z[i])&(V<=z[i]))/sum(U<=z[i]) + Remp[i]=sum((U>=1-z[i])&(V>=1-z[i]))/sum(U<=z[i]) + }
et pour la version théorique,
> Lg=(pcopula(copclayton,cbind(z,z)))/(z) > Rg=((1-2*(1-z)+pcopula(copclayton,cbind(1-z,1-z))))/(z) > plot(c(1-z,z),c(Lg,Rg))
De plus, on a des fonctions similaires pour la dépendance au sens faible, avec le code suivant pour la version théorique,
> Lg=log(pcopula(cop,cbind(z,z)))/log(z) > Rg=log((1-2*(1-z)+pcopula(cop,cbind(1-z,1-z))))/log(z) > Lg=1/Lg*2-1 > Rg=1/Rg*2-1
et celui là pour la version empirique
> z=seq(0,.5,by=.001) > v <- lossalae > U=rank(v[,1])/(nrow(v)+1) > V=rank(v[,2])/(nrow(v)+1) > Lemp=rep(NA,length(z)) > Remp=rep(NA,length(z)) > for(i in 1:length(z)){ + Lemp[i]=log(mean((U<=z[i])&(V<=z[i])))/log(mean(U<=z[i])) + Remp[i]=log(mean((U>=1-z[i])&(V>=1-z[i])))/log(mean(U<=z[i])) + } > Lemp=1/Lemp*2-1 > Remp=1/Remp*2-1
Bref, on peut utiliser ces fonctions sur des vrais échantillons. Considérons l’exemple classique loss-alae (où l’on couple les frais dans des sinistres assurés, et les frais payés par l’assureur). On souhaite ajuster une copule, sans trop savoir laquelle. On peut commencer par étudier la dépendance forte, et comparer avec une copule Gaussienne. La copule Gaussienne de référence possède ici le même rho de Spearman que l’échantillon dont on dispose,
> cor(lossalae,method="spearman") Loss ALAE Loss 1.000000 0.451872 ALAE 0.451872 1.000000 > library(copula) > paramgauss=.47 > paramclayton=.9 > paramgumbel=1.45 > copgauss=normalCopula(paramgauss) > copclayton=claytonCopula(paramclayton, dim = 2) > copgumbel=gumbelCopula(paramgumbel, dim = 2)
On obtient ici
La courbe verte est l’intervalle de confiance (ponctuel) à 95% pour une copule Gaussienne et un échantillon de même taille. On voit qu’on modélise mal la structure de dépendance. Avec une copule duale de Clayton, on obtient
et enfin pour une copule de Gumbel,
Bref, la copule de Gumbel semble réellement bien adaptée… Si on creuse en étudiant la dépendance au sens faible, on peut valider là aussi ce modèle. En effet, si la référence est la copule Gaussienne,
ou pour une copule de Clayton,
alors qu’une copule de Gumbel donnerait
