# Freakonometrics

## Tag - VaR

Thursday, September 27 2012

## Bounding sums of random variables, part 2

It is possible to go further, much more actually, on bounding sums of random variables (mentioned in the previous post). For instance, if everything has been defined, in that previous post, on distributions on $\mathbb{R}^+$, it is possible to extend bounds of distributions on $\mathbb{R}$. Especially if we deal with quantiles. Everything we've seen remain valid. Consider for instance two $\mathcal{N}(0,1)$ distributions. Using the previous code, it is possible to compute bounds for the quantiles of the sum of two Gaussian variates. And one has to remember that those bounds are sharp.

> Finv=function(u) qnorm(u,0,1)
> Ginv=function(u) qnorm(u,0,1)
> n=1000
> Qinf=Qsup=rep(NA,n-1)
> for(i in 1:(n-1)){
+ J=0:i
+ Qinf[i]=max(Finv(J/n)+Ginv((i-J)/n))
+ J=(i-1):(n-1)
+ Qsup[i]=min(Finv((J+1)/n)+Ginv((i-1-J+n)/n))
+ }

Actually, it is possible to compare here with two simple cases: the independent case, where the sum has a $\mathcal{N}(0,2)$ distribution, and the comonotonic case where the sum has a $\mathcal{N}(0,2^2)$ distribution.

>  lines(x,qnorm(x,sd=sqrt(2)),col="blue",lty=2)
>  lines(x,qnorm(x,sd=2),col="blue",lwd=2)

On the graph below, the comonotonic case (usually considered as the worst case scenario) is the plain blue line (with here an animation to illustrate the convergence of the numerical algorithm)

Below that (strong) blue line, then risks are sub-additive for the Value-at-Risk, i.e.

$H^{-1}(u)\leq F^{-1}(u)+G^{-1}(u)$

but above, risks are super-additive for the Value-at-RIsk. i.e.

$H^{-1}(u)\geq F^{-1}(u)+G^{-1}(u)$

(since for comonotonic variates, the quantile of the sum is the sum of quantiles). It is possible to visualize those two cases above, in green the area where risks are super-additive, while the yellow area is where risks are sub-additive.

Recall that with a Gaussian random vector, with correlation $r$ then the quantile is the quantile of a random variable centered, with variance $2(1+r)$. Thus, on the graph below, we can visualize case that can be obtained with this Gaussian copula. Here the yellow area can be obtained with a Gaussian copulas, the upper and the lower bounds being respectively the comonotonic and the countermononic cases.

But the green area can also be obtained when we sum two Gaussian variables ! We just have to go outside the Gaussian world, and consider another copula.

Another point is that, in the previous post, $C^-$ was the lower Fréchet-Hoeffding bound on the set of copulas. But all the previous results remain valid if $C^-$ is a lower bound on the set of copulas of interest. Especially

$\tau_{C^-,L}(F,G)\leq \sigma_{C,L}(F,G)\leq\rho_{C^-,L}(F,G)$

for all $C$ such that $C\geq C^-$. For instance, if we assume that the copula should have positive dependence, i.e. $C\geq C^\perp$, then

$\tau_{C^\perp,L}(F,G)\leq \sigma_{C,L}(F,G)\leq\rho_{C^\perp,L}(F,G)$

Which means we should have sharper bounds. Numerically, it is possible to compute those sharper bounds for quantiles. The lower bound becomes

$\sup_{u\in[x,1]}\left\{F^{-1}(u)+G^{-1}\left(\frac{x}{u}\right)\right\}$

while the upper bound is

$\sup_{u\in[0,x]}\left\{F^{-1}(u)+G^{-1}\left(\frac{x-u}{1-u}\right)\right\}$

Again, one can easily compute those quantities on a grid of the unit interval,

> Qinfind=Qsupind=rep(NA,n-1)
> for(i in 1:(n-1)){
+  J=1:(i)
+  Qinfind[i]=max(Finv(J/n)+Ginv((i-J)/n/(1-J/n)))
+  J=(i):(n-1)
+  Qsupind[i]=min(Finv(J/n)+Ginv(i/J))
+ }

We get the graph below (the blue area is here to illustrate how sharper those bounds get with the assumption that we do have positive dependence, this area been attained only with copulas exhibiting non-positive dependence)

For high quantiles, the upper bound is rather close to the one we had before, since worst case are probably obtained when we do have positive correlation. But it will strongly impact the lower bound. For instance, it becomes now impossible to have a negative quantile, when the probability exceeds 75% if we do have positive dependence...

> Qinfind[u==.75]
[1] 0

## Bounding sums of random variables, part 1

For the last course MAT8886 of this (long) winter session, on copulas (and extremes), we will discuss risk aggregation. The course will be mainly on the problem of bounding  the distribution (or some risk measure, say the Value-at-Risk) for two random variables with given marginal distribution. For instance, we have two Gaussian risks. What could be be worst-case scenario for the 99% quantile of the sum ? Note that I mention implications in terms of risk management, but of course, those questions are extremely important in terms of statistical inference, see e.g. Fan & Park (2006).

This problem, is sometimes related to some question asked by Kolmogorov almost one hundred years ago, as mentioned in Makarov (1981). One year after, Rüschendorf (1982) also suggested a proof of bounds calculation. Here, we focus in dimension 2. As usual, it is the simple case. But as mentioned recently, in Kreinovich & Ferson (2005), in dimension 3 (or higher), "computing the best-possible bounds for arbitrary n is an NP-hard (computationally intractable) problem". So let us focus on the case where we sum (only) two random variable (for those interested in higher dimension, Puccetti & Rüschendorf (2012) provided interesting results for a dual version of those optimal bounds).

Let $\Delta$ denote the set of univariate continuous distribution function, left-continuous, on $\mathbb{R}$. And $\Delta^+$ the set of distributions on $\mathbb{R}^+$. Thus, $F\in\Delta^+$ if $F\in\Delta$ and $F(0)=0$. Consider now two distributions $F,G\in\Delta^+$. In a very general setting, it is possible to consider operators on $\Delta^+\times \Delta^+$. Thus, let $T:[0,1]\times[0,1]\rightarrow[0,1]$ denote an operator, increasing in each component, thus that $T(1,1)=1$. And consider some function $L:\mathbb{R}^+\times\mathbb{R}^+\rightarrow\mathbb{R}^+$ assumed to be also increasing in each component (and continuous). For such functions $T$ and $L$, define the following (general) operator, $\tau_{T,L}(F,G)$ as

$\tau_{T,L}(F,G)(x)=\sup_{L(u,v)=x}\{T(F(u),G(v))\}$

One interesting case can be obtained when $T$is a copula, $C$. In that case,

$\tau_{C,L}(F,G):\Delta^+\times\Delta^+\rightarrow\Delta^+$

and further, it is possible to write

$\tau_{C,L}(F,G)(x)=\sup_{(u,v)\in L^{-1}(x)}\{C(F(u),G(v))\}$

It is also possible to consider other (general) operators, e.g. based on the sum

$\sigma_{C,L}(F,G)(x)=\int_{(u,v)\in L^{-1}(x)} dC(F(u),G(v))$

or on the minimum,

$\rho_{C,L}(F,G)(x)=\inf_{(u,v)\in L^{-1}(x)}\{C^\star(F(u),G(v))\}$

where $C^\star$ is the survival copula associated with $C$, i.e. $C^\star(u,v)=u+v-C(u,v)$. Note that those operators can be used to define distribution functions, i.e.

$\sigma_{C,L}(F,G):\Delta^+\times\Delta^+\rightarrow\Delta^+$

and similarly

$\rho_{C,L}(F,G):\Delta^+\times\Delta^+\rightarrow\Delta^+$

All that seems too theoretical ? An application can be the case of the sum, i.e. $L(x,y)=x+y$, in that case $\sigma_{C,+}(F,G)$ is the distribution of sum of two random variables with marginal distributions $F$ and $G$, and copula $C$. Thus, $\sigma_{C^\perp,+}(F,G)$ is simply the convolution of two distributions,

$\sigma_{C^\perp,+}(F,G)(x)=\int_{u+v=x} dC^\perp(F(u),G(v))$

The important result (that can be found in Chapter 7, in Schweizer and Sklar (1983)) is that given an operator $L$, then, for any copula $C$, one can find a lower bound for $\sigma_{C,L}(F,G)$

$\tau_{C^-,L}(F,G)\leq \tau_{C,L}(F,G)\leq\sigma_{C,L}(F,G)$

as well as an upper bound

$\sigma_{C,L}(F,G)\leq \rho_{C,L}(F,G)\leq\rho_{C^-,L}(F,G)$

Those inequalities come from the fact that for all copula $C$, $C\geq C^-$, where $C^-$ is a copula. Since this function is not copula in higher dimension, one can easily imagine that get those bounds in higher dimension will be much more complicated...

In the case of the sum of two random variables, with marginal distributions $F$ and $G$, bounds for the distribution of the sum $H(x)=\mathbb{P}(X+Y\leq x)$, where $X\sim F$ and $Y\sim G$, can be written

$H^-(x)=\tau_{C^- ,+}(F,G)(x)=\sup_{u+v=x}\{ \max\{F(u)+G(v)-1,0\} \}$

for the lower bound, and

$H^+(x)=\rho_{C^- ,+}(F,G)(x)=\inf_{u+v=x}\{ \min\{F(u)+G(v),1\} \}$

for the upper bound. And those bounds are sharp, in the sense that, for all $t\in(0,1)$, there is a copula $C_t$ such that

$\tau_{C_t,+}(F,G)(x)=\tau_{C^- ,+}(F,G)(x)=t$

and there is (another) copula $C_t$ such that

$\sigma_{C_t,+}(F,G)(x)=\tau_{C^- ,+}(F,G)(x)=t$

Thus, using those results, it is possible to bound cumulative distribution function. But actually, all that can be done also on quantiles (see Frank, Nelsen & Schweizer (1987)). For all $F\in\Delta^+$ let $F^{-1}$ denotes its generalized inverse, left continuous, and let $\nabla^+$ denote the set of those quantile functions. Define then the dual versions of our operators,

$\tau^{-1}_{T,L}(F^{-1},G^{-1})(x)=\inf_{(u,v)\in T^{-1}(x)}\{L(F^{-1}(u),G^{-1}(v))\}$

and

$\rho^{-1}_{T,L}(F^{-1},G^{-1})(x)=\sup_{(u,v)\in T^\star^{-1}(x)}\{L(F^{-1}(u),G^{-1}(v))\}$

Those definitions are really dual versions of the previous ones, in the sense that $\tau^{-1}_{T,L}(F^{-1},G^{-1})=[\tau_{T,L}(F,G)]^{-1}$ and $\rho^{-1}_{T,L}(F^{-1},G^{-1})=[\rho_{T,L}(F,G)]^{-1}$.

Note that if we focus on sums of bivariate distributions, the lower bound for the quantile of the sum is

$\tau^{-1}_{C^{-},+}(F^{-1},G^{-1})(x)=\inf_{\max\{u+v-1,0\}=x}\{F^{-1}(u)+G^{-1}(v)\}$

while the upper bound is

$\rho^{-1}_{C^{-},+}(F^{-1},G^{-1})(x)=\sup_{\min\{u+v,1\}=x}\{F^{-1}(u)+G^{-1}(v)\}$

A great thing is that it should not be too difficult to compute numerically those quantities. Perhaps a little bit more for cumulative distribution functions, since they are not defined on a bounded support. But still, if the goal is to plot those bounds on $[0,10]$, for instance. The code is the following, for the sum of two lognormal distributions $LN(0,1)$

> F=function(x) plnorm(x,0,1)
> G=function(x) plnorm(x,0,1)
> n=100
> X=seq(0,10,by=.05)
> Hinf=Hsup=rep(NA,length(X))
> for(i in 1:length(X)){
+ x=X[i]
+ U=seq(0,x,by=1/n); V=x-U
+ Hinf[i]=max(pmax(F(U)+G(V)-1,0))
+ Hsup[i]=min(pmin(F(U)+G(V),1))}

If we plot those bounds, we obtain

> plot(X,Hinf,ylim=c(0,1),type="s",col="red")
> lines(X,Hsup,type="s",col="red")

But somehow, it is even more simple to work with quantiles since they are defined on a finite support. Quantiles are here

> Finv=function(u) qlnorm(u,0,1)
> Ginv=function(u) qlnorm(u,0,1)

The idea will be to consider a discretized version of the unit interval as discussed in Williamson (1989), in a much more general setting. Again the idea is to compute, for instance

$\sup_{u\in[0,x]}\{F^{-1}(u)+G^{-1}(x-u)\}$

The idea is to consider $x=i/n$ and $u=j/n$, and the bound for the quantile function at point $i/n$ is then

$\sup_{j\in\{0,1,\cdots,i\}}\left\{F^{-1}\left(\frac{j}{n}\right)+G^{-1}\left(\frac{i-j}{n}\right)\right\}$

The code to compute those bounds, for a given $n$ is here

> n=1000
> Qinf=Qsup=rep(NA,n-1)
> for(i in 1:(n-1)){
+ J=0:i
+ Qinf[i]=max(Finv(J/n)+Ginv((i-J)/n))
+ J=(i-1):(n-1)
+ Qsup[i]=min(Finv((J+1)/n)+Ginv((i-1-J+n)/n))
+ }

Here we have (several $n$s were considered, so that we can visualize the convergence of that numerical algorithm),

Here, we have a simple code to visualize bounds for quantiles for the sum of two risks. But it is possible to go further...

Friday, February 10 2012

## exchangeability, credit risk and risk measures

Exchangeability is an extremely concept, since (most of the time) analytical expressions can be derived. But it can also be used to observe some unexpected behaviors, that we will discuss later on with a more general setting. For instance, in a old post, I discussed connexions between correlation and risk measures (using simulations to illustrate, but in the context of exchangeable risk, calculations can be performed more accurately). Consider again the standard credit risk problem, where the quantity of interest is the number of defaults in a portfolio. Consider an homogeneous portfolio of exchangeable risk. The quantity of interest is here

or perhaps the quantile function of the sum (since the Value-at-Risk is the standard risk measure). We have seen yesterday that - given the latent factor - (either the company defaults, or not), so that

i.e. we can derive the (unconditional) distribution of the sum

so that the probability function of the sum is, assuming that

Thus,  the following code can be used to calculate the quantile function
> proba=function(s,a,m,n){
+ b=a/m-a
+ choose(n,s)*integrate(function(t){t^s*(1-t)^(n-s)*
+ dbeta(t,a,b)},lower=0,upper=1,subdivisions=1000,
+ stop.on.error =  FALSE)$value + } > QUANTILE=function(p=.99,a=2,m=.1,n=500){ + V=rep(NA,n+1) + for(i in 0:n){ + V[i+1]=proba(i,a,m,n)} + V=V/sum(V) + return(min(which(cumsum(V)>p))) } Now observe that since variates are exchangeable, it is possible to calculate explicitly correlations of defaults. Here i.e. Thus, the correlation between two default indicators is then Under the assumption that the latent factor is beta distributed we get Thus, as a function of the parameter of the beta distribution (we consider beta distributions with the same mean, i.e. the same margin distributions, so we have only one parameter left, with is simply the correlation of default indicators), it is possible to plot the quantile function, > PICTURE=function(P){ + A=seq(.01,2,by=.01) + VQ=matrix(NA,length(A),5) + for(i in 1:length(A)){ + VQ[i,1]=QUANTILE(a=A[i],p=.9,m=P) + VQ[i,2]=QUANTILE(a=A[i],p=.95,m=P) + VQ[i,3]=QUANTILE(a=A[i],p=.975,m=P) + VQ[i,4]=QUANTILE(a=A[i],p=.99,m=P) + VQ[i,5]=QUANTILE(a=A[i],p=.995,m=P) + } + plot(A,VQ[,5],type="s",col="red",ylim= + c(0,max(VQ)),xlab="",ylab="") + lines(A,VQ[,4],type="s",col="blue") + lines(A,VQ[,3],type="s",col="black") + lines(A,VQ[,2],type="s",col="blue",lty=2) + lines(A,VQ[,1],type="s",col="red",lty=2) + lines(A,rep(500*P,length(A)),col="grey") + legend(3,max(VQ),c("quantile 99.5%","quantile 99%", + "quantile 97.5%","quantile 95%","quantile 90%","mean"), + col=c("red","blue","black", +"blue","red","grey"), + lty=c(1,1,1,2,2,1),border=n) +} e.g. with a (marginal) default probability of 15%, > PICTURE(.15) On this graph, we observe that the stronger the correlation (the more on the left), the higher the quantile... Note that the same graph can be plotted with on the X-axis the correlation, Which is quite intuitive, somehow. But if the marginal probability of default decreases, increasing the correlation might decrease the risk (i.e. the quantile function), > PICTURE(.05) (with the modified code to visualize the quantile as a function of the underlying default correlation) or even worse, > PICTURE(.0075) And it because all the more counterintuitive that the default probability decreases ! So in the case of a portfolio of non-very risky bond issuers (with high ratings), assuming a very strong correlation will lower risk based capital ! Thursday, February 9 2012 ## MAT8886 from tail estimation to risk measure(s) estimation This week, we conclude the part on extremes with an application of extreme value theory to risk measures. We have seen last week that, if we assume that above a threshold , a Generalized Pareto Distribution will fit nicely, then we can use it to derive an estimator of the quantile function (for percentages such that the quantile is larger than the threshold) It the threshold is , i.e. we keep the largest observations to fit a GPD, then this estimator can be written The code we wrote last week was the following (here based on log-returns of the SP500 index, and we focus on large losses, i.e. large values of the opposite of log returns, plotted below) > library(tseries) > X=get.hist.quote("^GSPC") > T=time(X) > D=as.POSIXlt(T) > Y=X$Close
> R=diff(log(Y))
> D=D[-1]
> X=-R
> plot(D,X)
> library(evir)
> GPD=gpd(X,quantile(X,.975))
> xi=GPD$par.ests[1] > beta=GPD$par.ests[2]
> u=GPD$threshold > QpGPD=function(p){ + u+beta/xi*((100/2.5*(1-p))^(-xi)-1) + } > QpGPD(1-1/250) 97.5% 0.04557386 > QpGPD(1-1/2500) 97.5% 0.08925095  This is similar with the following outputs, with the return period of a yearly event (one observation out of 250 trading days) > gpd.q(tailplot(gpd(X,quantile(X,.975))), 1-1/250, ci.type = + "likelihood", ci.p = 0.95,like.num = 50) Lower CI Estimate Upper CI 0.04172534 0.04557655 0.05086785  or the decennial one > gpd.q(tailplot(gpd(X,quantile(X,.975))), 1-1/2500, ci.type = + "likelihood", ci.p = 0.95,like.num = 50) Lower CI Estimate Upper CI 0.07165395 0.08925558 0.13636620  Note that it is also possible to derive an estimator for another population risk measure (the quantile is simply the so-called Value-at-Risk), the expected shortfall (or Tail Value-at-Risk), i.e. The idea is to write that expression so that we recognize the mean excess function (discussed earlier). Thus, assuming again that above (and therefore above that high quantile) a GPD will fit, we can write or equivalently If we substitute estimators to unknown quantities on that expression, we get The code is here > EpGPD=function(p){ + u-beta/xi+beta/xi/(1-xi)*(100/2.5*(1-p))^(-xi) + } > EpGPD(1-1/250) 97.5% 0.06426508 > EpGPD(1-1/2500) 97.5% 0.1215077  An alternative is to use Hill's approach (used to derive Hill's estimator). Assume here that , where is a slowly varying function. Then, for all , Since is a slowly varying function, it seem natural to assume that this ratio is almost 1 (which is true asymptotically). Thus i.e. if we invert that function, we derive an estimator for the quantile function which can also be written (which is close to the relation we derived using a GPD model). Here the code is > k=trunc(length(X)*.025) > Xs=rev(sort(as.numeric(X))) > xiHill=mean(log(Xs[1:k]))-log(Xs[k+1]) > u=Xs[k+1] > QpHill=function(p){ + u+u*((100/2.5*(1-p))^(-xiHill)-1) + } with the following Hill plot For yearly and decennial events, we have here > QpHill(1-1/250) [1] 0.04580548 > QpHill(1-1/2500) [1] 0.1010204 Those quantities seem consistent since they are quite close, but they are different compared with empirical quantiles, > quantile(X,1-1/250) 99.6% 0.04743929 > quantile(X,1-1/2500) 99.96% 0.09054039  Note that it is also possible to use some functions to derive estimators of those quantities, > riskmeasures(gpd(X,quantile(X,.975)),1-1/250) p quantile sfall [1,] 0.996 0.04557655 0.06426859 > riskmeasures(gpd(X,quantile(X,.975)),1-1/2500) p quantile sfall [1,] 0.9996 0.08925558 0.1215137 (in this application, we have assumed that log-returns were independent and identically distributed... which might be a rather strong assumption). An to conclude that post (and the session on extremes), here is the schedule for next weeks, Thursday, February 2 2012 ## quantile estimation (for rare events) During the course, we will discuss quantile estimation. We will use some financial series, > library(tseries) > X=get.hist.quote("^GSPC") > T=time(X) > D=as.POSIXlt(T) > Y=X$Close
> R=diff(log(Y))
> D=D[-1]
> plot(D,R)

In case it might not be possible to install some packages in the class room, it is also possible to use a local copy of the series,

>  base=read.table(
+  "http://freakonometrics.blog.free.fr/public/data/sp500.txt",
>  R=base$R > D=(as.Date(base$D))
>  plot(D,R)

Several libraries can be used to estimate standard risk measures, e.g. quantiles or expected shortfalls.

>  library(fExtremes)
>   VaR(R, alpha = 0.05, type = "sample", tail =  "upper")
95%
0.01725457
>  VaR(R, alpha = 0.005, type = "sample", tail =  "upper")
99.5%
0.04088206 

Or

> library(evir)
> riskmeasures(gpd(R,.015),.95)
p   quantile      sfall
[1,] 0.95 0.01743481 0.02742026
> riskmeasures(gpd(R,.015),.995)
p   quantile      sfall
[1,] 0.995 0.04079682 0.05598265

using that function, note that it is also possible to draw the log-likelihood function to visualize confidence intervals,

>  gpd.q(tailplot(gpd(R,.015)), .995, ci.type =
+ "likelihood", ci.p = 0.95,like.num = 50)
Lower CI   Estimate   Upper CI
0.03768945 0.04079682 0.04511344 

This week, we will see in details where the estimators of those quantities come from.

Friday, December 3 2010

## Millenium bridge, endogeneity and risk management

In less than 48 hours, last week two friends mentioned the Millennium Bridge as an illustration of a risk management concept. There are several documents with that example, here (for the initial idea of using the Millennium Bridge to illustrate issues in risk management) here or there, e.g.

When we mention resonance effects on bridges, we usually thing of the Tacoma Narrows Bridge (where strong winds set the bridge oscillating) or the Basse-Chaîne Bridge (in France, which collapsed on April 16, 1850, when 478 French soldiers marched across it in lockstep). In the first case, there is nothing we can do about it, but for the second one, this is why soldiers are required to break step on bridges.

But for the Millennium bridge, a 'positive feedback' phenomenon (known as Synchronous Lateral Excitation in physics) has been observed: the natural sway motion of people walking caused small sideways oscillations in the bridge, which in turn caused people on the bridge to sway in step, increasing the amplitude of the oscillations and continually reinforcing the effect. That has been described in a nice paper in 2005 (here). In the initial paper by Jon Danielsson and Hyun Song Shin, they note that "what is the probability that a thousand people walking at random will end up walking exactly in step? It is tempting to say “close to zero”, or “negligible”. After all, if each person’s step is an independent event, then the probability of everyone walking in step would be the product of many small numbers - giving us a probability close to zero. Presumably, this is the reason why Arup - the bridge engineers - did not take this into account. However, this is exactly where endogenous risk comes in. What we must take into account is the way that people react to their environment. Pedestrians on the bridge react to how the bridge is moving. When the bridge moves under your feet, it is a natural reaction for people to adjust their stance to regain balance. But here is the catch. When the bridge moves, everyone adjusts his or her stance at the same time. This synchronized movement pushes the bridge that the people are standing on, and makes the bridge move even more. This, in turn, makes the people adjust their stance more drastically, and so on. In other words, the wobble of the bridge feeds on itself. When the bridge wobbles, everyone adjusts their stance, which sets off an even worse wobble, which makes the people adjust even more, and so on. So, the wobble will continue and get stronger even though the initial shock (say, a gust of wind) has long passed. It is an example of a force that is generated and amplified within the system. It is an endogenous response. It is very different from a shock that comes from a storm or an earthquake which are exogenous to the system."

And to go further, they point out that this event is rather similar to what is observed in financial markets (here) by quoting The Economist from October 12th 2000 "So-called value-at-risk models (VaR) blend science and art. They estimate how much a portfolio could lose in a single bad day. If that amount gets too large, the VAR model signals that the bank should sell. The trouble is that lots of banks have similar investments and similar VAR models. In periods when markets everywhere decline, the models can tell everybody to sell the same things at the same time, making market conditions much worse. In effect, they can, and often do, create a vicious feedback loop. "

Friday, July 17 2009

## Quantile de sommes, ou somme de quantiles ?

Je reprends ici un commentaire que j'ai longtemps entendu dans mon expérience d'actuaire dans une vie antérieure (et qu'on peut lire - entre les lignes le plus souvent - dans certaines réactions de risk managers). Les quantiles sont la base des mesures de risques, en particulier via la VaR en finance, mais quand j'étais actuaire, on calculait des 90/10 pour désigner des quantiles à 90% des pertes. Donc tout ça n'est pas spécialement nouveau. Je rappelle que la fonction quantile est simplement  où
Les quantiles sont calculés par type de risque (en gestion des risques), et forcément à la fin, on souhaite un quantile de la somme des risques, et il a longtemps été coutume de prendre la somme des quantiles. Je vais donc discuter un peu le résultat suivant,
• L'égalité n'est pas vraie dans le cas indépendant
Suite à une espère de déformation professionnelle (à cause de la variance), beaucoup de monde pense que
mais ce résultat est faux.... Un exemple est d'ailleurs donné ci-dessous dans le cas de variables exponentielles.
• L'égalité est vraie dans le cas comonotone...
En fait, le résultat précédant est vrai dans le cas comonotone. La comonotonie correspond au cas de dépendance parfaite positive, de corrélation maximale (si la corrélation existe), i.e.  avec  une fonction strictement croissante. Alors
• ... mais ce n'est pas un "worst case"
Seconde bizarrerie, la comonotonie n'est pas le pire des cas. La comonotonie est souvent envisagée comme un worst case scenario mais il n'en est rien (cette idée se retrouve dans les QIS 3 par exemple, ici).
Formellement, si  (la classe de Fréchet), alors
et
En passant en dimension 2 (cas où la borne inférieure est effectivement une vraie fonction de répartition), si et  désignent respectivement des version anticomonotone et comonotone de , alors on peut se demander, en considérant une mesure de risque quelconque  si
En fait, Andre Tchen a montré en 1980 un résultat qui montre que la comonotonie est effectivement un worst case, mais seulement dans certains cas très limités.
Si est une fonction supermodulaire, i.e.

pour tout  et . Alors dans ce cas, pour tout ,
(la preuve se trouve ici). Je renvoie aux papiers de Michel Denuit, de Jan Dhaene ou de Marc Goovarts (ici ou ) pour des applications en tarification sur deux têtes, par exemple, ou sur les primes stop-loss en réassurance
• Borner (numériquement) le quantile d'une somme
Commençons par un exemple, avec deux lois exponentielles, car c'est le plus simple, et surtout on peut calculer explicitement les bornes (dans le cas général, on se contentera de méthode numérique comme l'avait fait Khalil ici).
Si  et , où, pour être plus précis,
est valide pour tout  quelle que soit la structure de dépendance entre et , où
Aussi, en prenant les inverses (car ces fonctions sont strictement croissante), on obtient une inégalité en terme de quantile, ou de Value-at-Risk,
pour tout .
Rappelons que dans le cas de variables indépendantes  alors que dans le cas comonotone . La figure ci-dessous montre les valeurs possible pour  la fonction de répartition de la somme, avec le cas indépendant en bleu, et le cas comotonone en rouge.On peut aussi visualiser les valeurs possibles pour les quantiles de la somme,Dans un cas général, Ce problème a été traité par Makarov en 1981, et par Frank, Nelsen et Schweizer en 1997 (ici), tout en restant en dimension 2. En fait, ils ont travaillé sur les bornes de la fonction de répartition d'une somme mais c'est pareil, comme on l'a vu juste auparavant. Ils ont montré le résultat suivant: si  alors pour tout
désigne la copule anticomonotone (qui est une copule en dimension 2),
en posant , et
On peut retrouver des résultats plus fin en creusant dans la direction de l'arithmetic probability, par exemple avec la thèse de Williamson (ici, cette thèse est remarquable, et trop peu citée, malheureusement). Je ne parle que de la dimension 2, l'extension en dimension supérieure est délicate (je renvoie aux travaux de Paul Embrechts sur le sujet, ici ou , par exemple).

Tuesday, March 17 2009

## Calculs de SCR, Solvency Capital Requirements

Pour reprendre le contexte général, Solvency II (l'analogue de la directive CRD pour les banques*) repose sur 3 piliers,

1. définir des seuils quantitatifs de calcul des provisions techniques des fonds propres, seuils qui seront à terme réglementaires, à savoir le MCR (Minimum Capital Requirement, niveau minimum de fonds propres en-dessous duquel l'intervention de l'autorité de contrôle sera automatique) et le SCR (Solvency Capital Requirement, capital cible nécessaire pour absorber le choc provoqué par une sinistralité exceptionnelle),
2. fixer des normes qualitatives de suivi des risques en interne aux sociétés, et définir comment l'autorité de contrôle doit exercer ses pouvoirs de surveillance dans ce contexte. Notons qu'en principe, les autorités de contrôle auront la possibilité de réclamer à des sociétés "trop risquées" de détenir un capital plus élevé que le montant suggéré par le calcul du SCR, et pourra les forcer àréduire leur exposition aux risques,
3. définir un ensemble d'information que les autorités de contrôle jugeront nécessaires pour exercer leur pouvoir de surveillance.

Cette histoire de pilliers peut s'illustrer de la manière suivante

Sur le premier pilier, assureurs et réassureurs devront mesurer les risques, et devront s'assurer qu'ils détiennent suffisement de capital pour les couvrir. En pratique, le CEIOPS et la Commission Européenne ont retenu une probabilité de ruine de 0,5%. Les calculs de capital se font alors de deux manières, au choix,
1. utiliser une formule standard. La formule ainsi que la calibration des paramètres ont été abordé à l'aide des QIS.
2. utiliser un modèle interne. Là dessus, le CEIOPS étudie les modalités d'évaluation.
En avril 2007, QIS3 a été lancé, afin de proposer une formule standard pour le calcul des MCR et SCR, en étudiant la problématique spécifique des groupes. En particulier, on trouve dans les documents la formule suivante (pour un calcul de basic SCR)
Cette formule sort du QIS3, mais on trouve des choses analogues dans Sandström (2004), par exemple,
Avec une contrainte forte sur la forme du SCR, il obtient alors
D'où sort cette formule ? Certains ont tenté des éléments de réponse, par exemple
Ce résultat n'est malheureusement pas très probant car il n'est jamais rien évoqué sur la dépendance entre les composantes, ce qui est troublant. Sandstôrm écrit quelque chose de similaire, même si pour lui "normalité" est ici entendu dans un cadre multivarié.
Une explication peut être trouvée dans un papier de Dietmar Pfeiffer et Doreen Straßburger (ici) paru dans le Scandinavian Actuarial Journal (téléchargeable ici). Il cherche à expliquer comment calculer le SCR,
Il note, et c'est effectivement l'intuition que l'on avait, que dans un monde Gaussien (multivarié), cette formule marche, aussi bien pour un SCR basé sur la VaR que la TVaR. En particulier, ils citent un livre de Sven Koryciorz, correspondant à sa thèse de doctorat, intitulée "Sicherheitskapitalbestimmung und –allokation in der Schadenversicherung. Eine risikotheoretische Analyse auf der Basis des Value-at-Risk und des Conditional Value-at-Risk", publiée en 2004.
Sinon, pour aller un peu plus loin, on peut aussi noter, dans les rapports du CEIOPS des déclarations un peu troublantes, par exemple
Il est pourtant facile de montrer que ce n'est pas le cas (même si c'est effectivement ce que préconise la "formule standard"). Le graphique ci-dessous montre l'évolution de la VaR d'une somme de risques corrélés (échangeables) en fonction de la corrélation sous-jacente: sur cet exemple, les risques très très corrélés sont moins risqués que des risques moyennement corrélés.
(la loi sous-jacente est une copule de Student). En revanche pour la TVaR, sur le même exemple, la TVaR de la somme est effectivement une fonction croissante avec la corrélation,

(plus de compléments dans les slides de l'école d'été à Lyon l'été dernier, ici).

* Pour reprendre des éléments de la page de wikipedia (ici), la directive européenne CRD (Capital Requirements Directive, i.e. Fonds Propres Réglementaires) transpose dans le droit européen les recommandations des accords de Bâle II, visant à calculer les fonds propres exigés pour les établissements financiers (i.e. directives 2006/48/CEet 2006/49/CE) .

Friday, December 19 2008

## Econométrie de la finance: calcul de VaR "instantannée"

Quelques lignes de code pour ceux qui n'auraient pas eu le temps de les taper, lors du TD du début de semaine

library(tseries); x=get.hist.quote("^FCHI")
y0=diff(log(x[,4]),1,1); yv=as.numeric(y0)
moy.ew=function(x,r){
m=rep(NA,length(x))
for(i in 1:length(x)){
m[i]=weighted.mean(x[1:i], rev(r^(0:(i-1))))}
return(m)}
sd.ew=function(x,r,m){
sd=rep(NA,length(x))
for(i in 1:length(x)){
sd[i]=weighted.mean((x[1:i]-m[i])^2, rev(r^(0:(i-1))))}
return(sd)}
quant=function(x,alpha){
q=rep(NA,length(x))
for(i in 201:length(x)){
q[i]=quantile(x[(i-200):i],alpha)}
return(q)}
VaR=quant(yv,.01)
plot(yv,ylab="")
moy.ew.y=moy.ew(yv,.95)
sd.ew.y=sqrt(sd.ew(yv,.95,moy.ew.y))
lines(moy.ew.y+qnorm(.01)*sd.ew.y,col="red")
moy.ew.y=moy.ew(yv,.5)
sd.ew.y=sqrt(sd.ew(yv,.5,moy.ew.y))
lines(moy.ew.y+qnorm(.01)*sd.ew.y,col="green")
lines(VaR,col="blue",lwd=2)

pour le code sur la volatilité locale - EWMA proposée par RiskMetrics.

Sinon quelques lignes de code sur l'approche par extrêmes, même si - en pratique - je modifie les fonctions tailplot et gpd.q de library(evir), pour ne pas avoir de tracé de graphique. De plus, le graphique ci-dessous est fait en continu (pas tous les 100 jours comme ci-dessous, pour gagner un peu de temps). Ce sont des seuils à -0.5% et -1% qui ont été retenu pour l'ajustement d'une loi de Pareto

library(evir)
quant.EVT=function(x,alpha,seuil){
q=rep(NA,length(x))
for(i in 1:100){
t=400+i*40
out=gpd(-x[(t-400):t],seuil)
OUT=tailplot(out,plot=FALSE)
Q=gpd.q(OUT,1-alpha)
q[400+i*40]=-Q[2]}
return(q)}