Finding Quantile Given Weird Density Function

For one of my Statistical Finance homeworks, I had to solve this question.

Given a density function \(f(x)\), where

\[f(x) = \int_{-\infty}^{\infty} \frac{1}{Z} \frac{|x+1|}{(x^2+1)^2} = 1\]

Find the quantile \(y\) that satisfies the \(F^{-1}(y) = 0.95\) where \(F(x)\) is the cumulative ditribution function of \(f(x)\) and \(F^{-1}\) is the quantile function (the inverse of \(F\)).

Obviously that’d mean that \(F(x) = \int f(x) \: dx\) and I should start by integrating \(f(x)\) and then finding its inverse. However, upon visual inspection, one should see that \(f(x)\) is a bitch to integrate.

Yes it is possible to integrate it by splitting the ranges for the absolute function, then splitting \((x+1) = x + 1\), then \(\frac{x}{(x^2+1)^2}\) can be integrated by substitution of \(u = x^2+1\) and the other portion can be done trigonometrically. Now because my integration-fu is too bad, let me solve for this numerically. So here’s the trick to (ab)use R.

First let’s solve for \(Z\). Since \(f(x)\) is a density function, then by definition of density functions, \(\int_{-\infty}^\infty f(x) \: dx = 1\). Then,

\[Z = \int_{-\infty}^{\infty} \frac{|x+1|}{(x^2+1)^2} = 1.785398\]

Great! Now let’s move to R. Unfortunately, R does not have a function for numerical computation of quantiles for arbitrary distribution functions. However, we can build one.

Finding \(y\) above is equivalent to the following optimization problem:

\[y^* = \arg\min_{y} (F(y) - 0.95)^2\]

Then all we need is an optimization routine to find \(y\) that minimizes the squared error \(F(y) - 0.95)^2\). In R there exists the library nlminb for numerical optimization. The functions are below:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
CDF = function (x, dist) {
integrate(f = dist, lower = -Inf, upper = x)$value
}
objective = function (x, quantile, dist) {
(CDF(x, dist) - quantile)^2
}
find_quantile = function (dist, quantile) {
result = nlminb(start = 0, objective = objective, quantile = quantile, dist = dist)$par
return (result)
}
crazy_eqn = function (x) {
abs(x+1)/((x^2+1)^2)
}
z = integrate(crazy_eqn, lower=-Inf, upper=Inf)
crazy_eqn_to_integrate = function(x) {
z$value * (abs(x+1)/((x^2+1)^2))
}
y = find_quantile(dist = crazy_eqn_to_integrate, quantile = 0.95)

Using this, we find that \(y = 0.03161596\). Let’s check this!

1
2
integrate(crazy_eqn_to_integrate, lower=-Inf, upper=alpha)
0.95 with absolute error < 1.8e-11

Booyah. This is why I’m not a statistics major.