Taleb: Statistics and Risk Two Different Businesses

Towards the end of this section, Taleb inserts a sidebar as follows:

Consider the right tail $K^{+} \textup{ in } \mathbb{R}^{+}$ and the left tail $K^{-} \textup{ in } \mathbb{R}^{-}$ . Without specifying the support of the distribution:
Definition 1.3 (Probability swamps payoff (thin tails)).
$\lim_{K^{+ } \rightarrow \infty } E \left [ X \mid_{X > K^{+}} \right ]=K^{+} \textup{ and } \lim_{K^{- } \rightarrow -\infty } E \left [X \mid_{X < K^{-}} \right ]=K^{-}$
Definition 1.4 (Payoff swamps probability (fat tails)).
$\exists \textup{ } \lambda^{+} > 1 \textup{ or } \lambda^{-} > 1 \textup{ s.t.}$
$\lim_{K^{+ } \rightarrow \infty }E \left [ X \mid_{X > K^{+}} \right ]=\lambda^{+}. K^{+} \textup{ and } \lim_{K^{- } \rightarrow -\infty }E \left [ X \mid_{X < K^{-}} \right ]= \lambda^{-}.K^{-}$

N.B: I added the “+” and “-” superscripts to K in the last line where they were missing.

I was curious to look at this behavior for a set of distributions: normal, exponential, gamma, weibull, pareto, and cauchy. I created some functions which I provide below which will allow exploration of the effect of the parameters of the distributions. I also provide a script that can be used to run a set of examples, one for each distribution. The interesting ones in terms of the sidebar above are pareto and cauchy, but the differences between the others are of interest too.

I am not sure I understand what is special about the cases where the constant of proportionality, λ > 1. Do we say payoff swamps probability if λ = 1.1? or 5? It’s the transition from “Probability swamps Payoff” to “Payoff swamps Probability” that occurs if λ > 1 that I am not understanding. I am guessing that if a risk manager were to casually assume that his expected loss would approximate K, then he might be in for a shock if it is actually 2.K!

There’s a reference to Chapter 4 which I have yet to read – perhaps I will “get it” at that time.

Conditional Pay-Off

I am going to work with the following definition of pay-off throughout and work with skewed distributions with the tail out to the right. The same principles apply if the tail is out to the left and we talk about an expected value between – ∞ and κ.

$\displaystyle E \left [ X \mid_{x > \kappa} \right ]=\frac{\int_{\kappa }^{\infty }x.f_{X}(x)dx}{\int_{\kappa }^{\infty}f_{X}(x)dx}$

where, f_X(x) is the probability density function of random variable “X”. I am only going to look at continuous distributions.

Cauchy Distribution

pdf:

$\displaystyle f_{X}(x; \mu, \sigma)= \frac{1}{\pi \sigma \left ( 1 + \left ( \frac{x-\mu}{\sigma} \right )^{2} \right )}$

Wikipedia says the following:

The Cauchy distribution is often used in statistics as the canonical example of a “pathological” distribution since both its mean and its variance are undefined.

So now I am curious. How does it behave in the tail? Let’s see …

$\displaystyle E=E \left [ X \mid x > \kappa;0,\sigma \right ]=\frac{\int_{\kappa}^{\infty} \frac{x}{\pi \sigma(1+(x/ \sigma)^{2})}dx}{\int_{\kappa}^{\infty} \frac{1}{\pi \sigma(1+(x/ \sigma)^{2})}dx}= \frac{N}{D}$

We need only investigate the numerator …

$\displaystyle N=\int_{\kappa}^{\infty}\frac{x}{\pi\sigma(1+(x/\sigma)^{2})}dx=\left [ \frac{\sigma \ln (1+x^{2})}{2\pi} \right ]_{\kappa}^{\infty}=\infty$

So E is infinite, and so is E / κ. So I guess λ is infinite?

Normal Distribution

pdf:

$\displaystyle f_{X}(x;\mu,\sigma^{2})=\frac{1}{\sigma.\sqrt{2.\pi}}\exp(-((x-\mu)^{2}/2.\sigma^{2}))$

Expected pay-off:

$\displaystyle E = E \left [ X \mid x > \kappa;\mu, \sigma \right ]=\mu + E \left [X \mid x > k;0, \sigma \right ], k=\kappa - \mu$

Translating the function so the mean is zero, and adding it back later.

$\displaystyle E[X\mid x > k;0,\sigma^{2}]=\frac{\int_{k}^{\infty}\frac{x.\exp(-x^{2}/2.\sigma^{2})}{\sigma.\sqrt{2.pi}}dx}{1 - F(k;0,\sigma^{2})} = \frac{N}{D}$

Where F(k) is the cdf of the distribution, so 1 – F(k) is the area under the right tail of the pdf.
Integrating the numerator …

$\displaystyle N=\frac{-\sigma}{\sqrt{2.\pi}}\left [ \exp\left (\frac{-x^{2}}{2.\sigma^{2}}\right ) \right ]_{k}^{\infty}=\frac{\sigma.\exp(-k^{2}/2.\sigma^{2})}{\sqrt{2.\pi}}$

Finally:

$\displaystyle E=\mu+\frac{\sigma.\exp(-k^{2}/2.\sigma^{2})}{\sqrt{2.\pi}.(1 - F(k;0,\sigma^{2}))}=\mu+\frac{\sigma.\exp(-(\kappa - \mu)^{2}/2.\sigma^{2})}{\sqrt{2.pi}.(1 - F(\kappa;\mu,\sigma^{2}))}$

Now, to illustrate Taleb’s point we can divide E / κ and see what happens as κ gets larger and larger:

E / κ = λ is clearly heading to one. Since the mean is 3 and sd is 2, the value on the x-axis of 12 is equivalent to 4.5 standard deviations from the mean. This gives a sense of how quickly the normal distribution approaches the condition λ = 1. It very rapidly approaches 1 + μ / κ. If μ were zero, the approach to λ = 1 would be even quicker.

Exponential Distribution

pdf:

$\displaystyle f_{X}(x;\lambda)=\lambda.\exp(-\lambda.x)$

Expected pay-off:

$\displaystyle E \left [ X \mid x > \kappa;\lambda \right ]=\frac{\int_{\kappa}^{\infty}x.\lambda.\exp(-\lambda.x)dx}{\int_{\kappa}^{\infty}\lambda.\exp(-\lambda.x)dx}=\frac{N}{D}$

Integrating by parts:
$\displaystyle v=x, du=\lambda.e^{-\lambda.x}.dx \Rightarrow N=\left [ -x.e^{-\lambda.x}\right ]_{\kappa}^{\infty}-\int_{\kappa}^{\infty}-e^{-\lambda.x}dx=\left ( \kappa+\lambda \right ).e^{-\lambda.\kappa}$

and

$\displaystyle D=\left [ -e^{-\lambda.x}\right ]_{\kappa}^{\infty}=e^{-\lambda.\kappa}$

$\displaystyle E=E\left [ X\mid x> \kappa;\lambda \right ]=\frac{N}{D}=\kappa+\lambda$

The chart at the right illustrates how E / κ = 1 + λ / κ behaves with κ.

Obviously, as κ increases E / κ tends to 1. (Note the “λ” here is NOT Taleb’s λ it is the rate constant in the exponential distribution). In this case the rate, λ, was 2, so the standard deviation would be 1/2. As you would expect, E / κ is approaching 1 more slowly than for the normal distribution which contains the exp(-x^2) term rather than exp(-x). A value of κ = 7 is 13 standard deviations above the mean (1/2).

Gamma Distribution

pdf:

$\displaystyle f_{X}(x;\alpha,\beta)=\frac{\beta^{\alpha}.x^{\alpha-1}.e^{-\beta.x}}{\Gamma(\alpha)}; \alpha=shape,\beta=rate$

Expected pay-off:

$\displaystyle E\left [ X\mid x> \kappa \right ]=\frac{\int_{\kappa}^{\infty}\beta^{\alpha}.x^{\alpha}.e^{-\beta.x}dx/\Gamma(\alpha)}{1-F(\kappa;\alpha,\beta)}=\frac{N}{D}$

Where F(k) is the cdf of the distribution, so 1 – F(k) is the area under the right tail of the pdf.

Integrating the numerator by parts:

$\displaystyle v=x^{\alpha}, du=\frac{\beta^{\alpha}.e^{-\beta.x}.dx}{\Gamma(\alpha)} \Rightarrow$

$\displaystyle N=\left [ \frac{-\beta^{\alpha-1}.x^{\alpha}.e^{-\beta.x}}{\Gamma(\alpha)} \right ]_{\kappa}^{\infty}-\frac{\alpha}{\beta}\int_{\kappa}^{\infty}\frac{-\beta^{\alpha}.x^{\alpha-1}.e^{-\beta.x}dx}{\Gamma(\alpha)}$

$\displaystyle = \frac{\beta^{\alpha-1}.\kappa^{\alpha}.e^{-\beta.\kappa}}{\Gamma(\alpha)} +\frac{\alpha.D}{\beta}$

so, substituting θ for 1/β so we use scale rather than rate …

$\displaystyle E=E \left [ X \mid x > \kappa;\alpha,\theta \right ]=\frac{\kappa^{\alpha}.e^{-\kappa/\theta}}{\Gamma(\alpha).\theta^{\alpha-1}.\left (1-F(\kappa;\alpha,\theta)\right )} +\alpha.\theta$

Once again we get the performance Taleb promised, with λ = E / κ approaching 1 as κ increases. With the parameters in this example the mean should be 3.75 (scale.shape) and the standard deviation 2.4 (scale.√(shape)). So, 10 standard deviations is approaching 30 on the x-axis. This distribution is converging on 1 at a slower rate than the normal. Given that the gamma probability distribution contains a x^α.e^(-x) term compared to the exponential distribution’s e^(-x) term, I would expect E / κ to approach 1 even more slowly than exponential. I would need to figure out a basis on which to compare this aspect of each distribution (probably compare limit as κ -> ∞ of d(E / κ)/dκ for each distribution).

Weibull Distribution

pdf:

$\displaystyle f_{X}(x;\lambda,k)=\frac{k}{\lambda}\left (\frac{x}{\lambda} \right )^{k-1}\exp(-(x/\lambda)^{k}); x \geq 0, k=shape,\lambda=scale$

Expected pay-off:

$\displaystyle E=E\left [ X \mid x > \kappa;\lambda,k \right ]=\frac{\int_{\kappa}^{\infty}\frac{k}{\lambda}\left (\frac{x}{\lambda} \right )^{k-1} x. \exp(-(x/\lambda)^{k})dx}{1-F(\kappa;\lambda,k)}=\frac{N}{D}$

Substitute

$\displaystyle t=(x/\lambda)^{k}$ :

$\displaystyle x=\lambda t^{1/k}, dx=\frac{\lambda t^{1/k - 1}dt}{k} \Rightarrow$

$\displaystyle N=\lambda\int_{\kappa}^{\infty}t^{1/k}e^{-t}dt=\lambda\Gamma(1+1/k,\kappa )$

Where Γ(α, β) is the upper incomplete gamma function. In “R” one can access the incomplete gamma functions via the gamma cumulative distribution function (pgamma(x, shape, scale, lower.tail)):

$\displaystyle 1-F(x;kk, \theta)=\frac{\Gamma(kk, x/\theta)}{\Gamma(kk)}, \theta=1 \Rightarrow$
$\displaystyle \Gamma(kk, x)=\Gamma(kk)(1-F(x;kk, 1))$

Substituting:

$\displaystyle kk=1+1/k \Rightarrow E=\frac{\Gamma(1+1/k)(1-F(\kappa;1+1/k, 1))}{1-F(\kappa;\lambda,k)}$

The only problem is this doesn’t seem to work – if anyone can see what I have wrong please let me know.

In the mean time I created a numerical integration approach which seems to work but is not going to be accurate. It generates the adjacent chart for E / κ. It appears to be approaching 1 in an asymptotic fashion, but I can’t be sure until I get the equations above working!

Pareto Distribution

pdf:

$\displaystyle f_X(x;\mu,\sigma,\gamma,\alpha)=$

$\displaystyle \frac{\alpha\sigma^{-1/\gamma}}{\gamma} \left ( x-\mu \right )^{\frac{1}{\gamma}-1} \left ( \left ( \frac{\sigma}{x-\mu} \right )^{-1/\gamma}+1\right )^{\left ( -\alpha-1 \right )}; x> \mu$

Expected pay-off:
$\displaystyle E=E\left [X \mid x > \kappa; \mu,\sigma,\gamma,\alpha \right ] =$

$\displaystyle \frac{\int_{\kappa}^{\infty} \frac{x \alpha \sigma^{-1/\gamma}}{\gamma} \left ( x-\mu \right )^{\frac{1}{\gamma}-1} \left ( \left ( \frac{\sigma}{x-\mu} \right )^{-1/\gamma}+1\right )^{\left ( -\alpha-1 \right )}dx}{1-F(\kappa; \mu,\sigma,\gamma,\alpha)}$

$\displaystyle =\frac{N}{D}$

Since I haven’t figured out that integration yet, as with my previous posts, I resorted to a numerical integration with the results of E / κ in the chart at the right. I cannot say for sure as I need an analytic solution, but it sure looks like E / κ is heading for around 2 rather than 1! Since I used inequality, γ = 0.75 and tail index, α = 1.5, i.e. α<= 2.γ => the variance of the distribution is undefined (fat-tailed). So we have an example of a fat tailed distribution that confirms Taleb’s point.

If we raise the inequality, gamma=0.9, as in the second chart, E / kappa seems to be approaching a limit higher than with gamma = 0.75.

Going the other way, making inequality, gamma=0.5, since variance is now defined (alpha > 2.gamma) we ought to now have a thin-tailed distribution and E / kappa should tend towards 1. But it doesn’t appear to – see the chart to the right. I conclude that I must not yet understand Taleb’s definition of fat-tailed vs thin-tailed. More on this later once I figure it out!

Here is the R script containing all the functions needed to explore the material above:

lambdaNormal &lt;- function(mean=0, sd=1){
	xSeq &lt;- seq(qnorm(.001, mean, sd), qnorm(0.999, mean, sd), 0.001)
	dev.new()
	plot(xSeq, dnorm(xSeq, mean, sd), t="l", main=paste0("Normal pdf (mean=", mean, " st. dev.=",
         sd, ")"), xlab="X", ylab="Probability density")
	
	kSeq &lt;- seq(qnorm(.75, mean, sd), qnorm(0.999999, mean, sd), 0.01)
	eNorm &lt;- mean + sd * exp(-(kSeq - mean)^2/2/sd^2) / sqrt(2*pi) / pnorm(kSeq, mean, sd,
         lower.tail=F)
	normLambda &lt;- eNorm / kSeq
	dev.new()
	plot(kSeq, normLambda, t="l", main="Expected Tail Value (x &gt;= k) / Tail (k)",
         sub=paste0("Distribution: Normal (mean=", mean, ", std dev=", sd, " )"),
         ylab="E[X|x&gt;k] / k", xlab="k", xlim=c(min(kSeq), max(kSeq)),
         ylim=c(0.95, max(normLambda)))
	abline(h=1, col="red")
}
lambdaExponential &lt;- function(rate=1){
	xSeq &lt;- seq(qexp(.001, rate), qexp(0.999, rate), 0.001)
	dev.new()
	plot(xSeq, dexp(xSeq, rate), t="l", main=paste0("Exponential pdf (rate=", rate, " )"),
         xlab="X", ylab="Probability Density")
	kSeq &lt;- seq(qexp(.75, rate), qexp(0.999999, rate), 0.01)
	eExp &lt;- kSeq + rate
	expLambda &lt;- eExp / kSeq
	dev.new()
	plot(kSeq, expLambda, t="l", main="Expected Tail Value (x &gt;= k) / Tail (k)",
         sub=paste0("Distribution: Exponential (rate=", rate, " )"), ylab="E[X|x&gt;k] / k",
         xlab="k", xlim=c(min(kSeq), max(kSeq)), ylim=c(0.95, max(expLambda)))
	abline(h=1, col="red")
}
lambdaGamma &lt;- function(shape=1, scale=1){
	xSeq &lt;- seq(qgamma(.001, shape, scale=scale), qgamma(0.999, shape, scale=scale), 0.001)
	dev.new()
	plot(xSeq, dgamma(xSeq, shape=shape, scale=scale), t="l",
         main=paste0("Gamma pdf (shape=", shape, ", scale=", scale, " )"),
         xlab="X", ylab="Probability Density")
	kSeq &lt;- seq(qgamma(.75, shape, scale=scale), qgamma(0.999999, shape, scale=scale), 0.01)
	eGamma &lt;- shape * scale + (kSeq^shape) * exp(-kSeq / scale) / gamma(shape) /
         scale^(shape - 1) / pgamma(kSeq, shape, scale=scale, lower.tail=F)
	gammaLambda &lt;- eGamma / kSeq
	dev.new()
	plot(kSeq, gammaLambda, t="l", main="Expected Tail Value (x &gt;= k) / Tail (k)",
         sub=paste0("Distribution: Gamma (shape=", shape, ", scale=", scale, " )"),
         ylab="E[X|x&gt;k] / k", xlab="k", xlim=c(min(kSeq), max(kSeq)), ylim=c(0.95, max(gammaLambda)))
	abline(h=1, col="red")
}
eTailPareto4 &lt;- function(x, increment=0.001, mu=0, sigma=1, gamma=0.75, alpha=1.5){
	expectedValue &lt;- ePareto4(mu, sigma, gamma, alpha, T)
	xRange &lt;- seq(0, x, increment)
	ifelse(is.na(expectedValue), NA, (expectedValue - increment * sum(xRange * dPareto4(xRange, 
         mu, sigma, gamma, alpha, T))) / (1 - pPareto4(x, mu, sigma, gamma, alpha, T)))
}
lambdaPareto4 &lt;- function(mu=0, sigma=1, gamma=0.75, alpha=1.5){
# Pareto4 - have to do integration numerically
	xSeq &lt;- seq(mu-1, qPareto4(p=0.95, mu, sigma, gamma, alpha, T), 0.001)
	dev.new()
	plot(xSeq, dPareto4(xSeq, mu, sigma, gamma, alpha, T), t="l",
         main=paste0("Pareto pdf (loc.=", mu, ", scale=", sigma, ", inequality=",
         gamma, ", tail index=", alpha, " )"), xlab="X", ylab="Probability Density")
	if (gamma &lt; alpha){
		kSeq=seq(qPareto4(0.9, mu, sigma, gamma, alpha, T),
                 qPareto4(0.999, mu, sigma, gamma, alpha, T), 0.001)
		ePareto &lt;- sapply(kSeq, eTailPareto4, sigma=sigma, gamma=gamma, alpha=alpha)
		pareto4Lambda &lt;- ePareto / kSeq
		dev.new()
		plot(kSeq, pareto4Lambda, t="l", main="Expected Tail Value (x &gt;= k) / Tail (k)",
                 sub=paste0("Distribution: Pareto (loc.=", mu, ", scale=", sigma, ",
                 inequality=", gamma, ", tail index=", alpha, " )"),
                 ylab="E[X|x&gt;k] / k", xlab="k", xlim=c(min(kSeq), max(kSeq)),
                 ylim=c(min(pareto4Lambda) * 0.9, max(pareto4Lambda)))
	}else{
		message("Gamma is greater than or equal to alpha - loss is undefined")
	}
}
eTailWeibull &lt;- function(x, increment=0.001, shape=1, scale=1){
	expectedValue &lt;- scale * gamma(1 + 1 / shape)
	xRange &lt;- seq(0, x, increment)
	(expectedValue - increment * sum(xRange * dweibull(xRange, shape, scale))) /
         pweibull(x, shape, scale, lower.tail=F)
}
igamma &lt;- function(s, x=0, upper=T) gamma(s) * pgamma(x, shape=s, scale=1, lower.tail=!upper)
lambdaWeibull &lt;- function(shape=1, scale=1){
	xSeq &lt;- seq(0, qweibull(0.999, shape, scale), 0.001)
	dev.new()
	plot(xSeq, dweibull(xSeq, shape, scale), t="l",
         main=paste0("Weibull pdf (shape=", shape, ", scale=", scale, " )"),
         xlab="X", ylab="Probability Density")
	kSeq &lt;- seq(qweibull(0.99, shape, scale), qweibull(0.9999, shape, scale), 0.0001)
# This works ...evalue at 2 ~ 2.3. at 10 ~10.6
	neWeibull &lt;- sapply(kSeq, eTailWeibull, shape=shape, scale=scale)
	nWeibullLambda &lt;- neWeibull / kSeq
	dev.new()
	plot(kSeq, nWeibullLambda, t="l", main="Expected Tail Value (x &gt;= k) / Tail (k)",
         sub=paste0("Distribution: Numeric Weibull (shape=", shape, ", scale=", scale, " )"),
         ylab="E[X|x&gt;k] / k", xlab="k", xlim=c(min(kSeq), max(kSeq)),
         ylim=c(min(nWeibullLambda) * 0.9, max(nWeibullLambda)))
	abline(h=1, col="red")
	
}
# Show that as k -&gt; +infinity, E[X | x &gt;= K+] -&gt; Lambda.K+ for fat tails, but -&gt; k+ for thin tailed
# List which distributions to view: normal, exponential, gamma, pareto4, weibull
view &lt;- c()

if ("normal" %in% view){
# Normal
	mu &lt;- 0
	stdDev &lt;- 2
	lambdaNormal(mu, stdDev)
}
if("exponential" %in% view){
# Exponential dist
	lambda &lt;- 2
	lambdaExponential(lambda)
}
if ("gamma" %in% view){
# Gamma:
	shape &lt;- 5/2
	scale &lt;- 1.5
	lambdaGamma(shape, scale)
}
if("pareto4" %in% view){
# Pareto4 - have to do integration numerically
	location &lt;- 0
	scale &lt;- 1
	inequality &lt;- 0.75
	tailIndex &lt;- 1.5
	lambdaPareto4(location, scale, inequality, tailIndex)
}
if("weibull" %in% view){
# Weibull
	shape &lt;- 1.5
	scale &lt;- 2
	lambdaWeibull(shape, scale)
}