[SOLVED] CS python finance project

project

Course Project: An Arbitrage-Free Smile Interpolator
Objectives
Implement an arbitrage free smile interpolator SmileAF.
Use the arbitrage free smile interpolator to construct local volatility model
Use PDE with local volatility model to price a given set of European options (strike in delta $times$ maturity)
Compare the price errors of arbitrage-free smile interpolator and the cubic spline smile interpolator

Smile Arbitrage
European call prices are monotonically decreasing with respect to the strike:

begin{align}
C(S_0, K_1, T, sigma(K_1), r, q) geq C(S_0, K_2, T, sigma(K_2), r, q) ~text{for}~K_1 < K_2 end{align}The European call price as a function of strike has to be convex every where: for any three points $K_1 < K_2 < K_3$begin{align} frac{C(K_2) – C(K_1) } {K_2 – K_1} < frac{C(K_3) – C(K_2) } {K_3 – K_2} end{align}orbegin{align} C(K_2)< C(K_3)frac{K_2 – K_1} {K_3 – K_1} + C(K_1)frac{K_3-K_2} {K_3 – K_1} end{align}This is also equivalent to “butterfly price has to be non-negative”.When Could Smile Arbitrage Happen?The undiscounted call price is the expectation of payoff under risk neutral measurebegin{align}C(K) = E[max(S-K, 0)]end{align}And expectation is an integral over the probability density function $p(s)$begin{align}C(K) = int_{K}^{+infty} (s-K) p(s) dsend{align}The 1st non-arbitrage condition translates tobegin{align}& C(K_1) – C(K_2) = left[ int_{K_1}^{K_2} (s-K_1) p(s) ds+ int_{K_2}^{+infty} (K_2-K_1) p(s) ds right]end{align}which is positive by definition if $K_2 > K_1$.

The 2nd non-arbitrage condition translates to

begin{align}
& C(K_3)frac{K_2 K_1} {K_3 K_1} + C(K_1)frac{K_3-K_2} {K_3 K_1} C(K_2) \
= & frac{K_3-K_2} {K_3 K_1} int_{K_1}^{K_2} (s-K_1) p(s) ds + frac{K_2 K_1} {K_3 K_1} int_{K_2}^{K_3} (K_3 s) p(s) ds
end{align}which is also positive by definition if $K_3 > K_2 > K_1$.

So, when could smile arbitrage happen? When the probability density does not exist. If we can start with valid probability density function $p(s)$, arbitrage-freeness is guaranteed by construction.

Arbitrage Free Smile (Based on [Fengler 2009])
We consider smile construction for a given expiry $T$.
Start with $N$ discrete sample strike points

begin{align}
vec{k} = [k_1, k2, ldots, k{N}]^{top}
end{align}

Try to solve for undiscounted call prices for these $N$ sample points

begin{align}
vec{c} = [c_1, c_2, ldots, c_N]^{top}
end{align}

For the undiscounted call price $C(K)$ for any $K$, we can interpolate using cubic spline over the sample points $(k_i, c_i)$. (Note that we are using cubic spline to interpolate the prices, not volatility)

The second derivative of call price with respect to strike is the probability density function:

begin{align}
frac{d C}{d K} & =d frac{int_K^{infty} Sp(S) dS}{dK} d frac{Kint_K^{infty} p(S) dS}{dK} = -Kp(K) left( int_K^{infty} p(S) dS K p(K)right) = -int_K^{infty} p(S) dS
frac{d^2 C}{d K^2} & = p(K)
end{align}

So $c_i$ is probability density function at $k_i$, we denote it as $p_i$

Second derivatives in cubic spline interpolation form line segments. Cubic spline on $C(K)$, means linearly interpolate on probability density. If $p_i$ are all positive, the whole pdf is positive by construction no smile arbitrage.

For tails call prices are almost linear if strike is very far away from spot, we can use natural cubic spline: $p_1 = p_N = 0$.

Our problem is to solve for $[c_1, c_2, ldots, c_{N}, p_2, ldots, p_{N-1}]$

Inputs to our problem
Same as our Cubic Spline smile interpolator, we have the input marks to start with to construct the Arb-Free(AF) smile interpolator:

Marks:strike to volatility pairs, denote as $(hat k_j, sigma_j)$, for $j in [1, 2, ldots, M]$. In our case, $M=5$.

We would like to match the marks exactly. And we cannot directly construct a cubic spline using the $M$ points too coarse and distribution is not realistic.

Problem Definition
We use $N = 50$ sample points, ranging from $[k_1 = S e^{(r_d r_f)T -frac12sigma_{ATM}^2T 5 sigma_{ATM} sqrt{T}}, k_N = S e^{(r_d r_f)T -frac12sigma_{ATM}^2T + 5 sigma_{ATM} sqrt{T} }]$, i.e., $pm 5$ standard deviation based on $sigma_{ATM}$.
$sigma_{ATM}$ is implied volatility of the middle point of the input marks.
We also assume the strike of the middle point of the input marks is the forward ATM forward convention.
The sample points are equally spaced, denote the length of the segment $u = frac{k_N k_1}{N-1}$
We would like the call prices to be as smooth as possible minimize the change of the slopes
We want to match exactly the $M$ input marks.
This is a constrained optimization problem.

Constraints

Cubic spline interpolation imposes the constraints that the left and right first derivative of a point have to match, it can be derived by matching the firstderivative of the left and right segments for point $i$ we have the condition

begin{align}
c{i+1} + c{i-1} 2 c_{i} = (frac23 pi + frac16 p{i+1} + frac16 p_{i-1}) u^2
end{align}

The cubic spline constraints translate to the linear system

begin{align} underbrace{begin{pmatrix}
1 & -2 & 1 & 0 & ldots & 0
0 & 1 & -2 & 1 & ddots & vdots
vdots & ddots & ddots & ddots & ddots & 0
0 & ldots & 0 & 1 & -2 & 1
end{pmatrix}}{vec{Q}{(N-2) times N}}
begin{pmatrix}
c_1
c_2
vdots
cN
end{pmatrix} =
underbrace{u^2
begin{pmatrix}
frac23 & frac16 & 0 & ldots & 0
frac16 & frac23 & frac16 & ddots & vdots
0 & ddots & ddots & ddots & 0
vdots & ddots & frac 1 6 & frac23 & frac16
0 & ldots & 0 & frac 1 6 &frac23
end{pmatrix}}{vec{R}_{(N-2) times (N-2)}}
begin{pmatrix}
p_2
p3
vdots
p{N-1}
end{pmatrix}
end{align}

If we define
begin{align}
vec{x} =
begin{pmatrix}
vec{c}^{top}
vec{p}^{top}
end{pmatrix}, ~~~
vec{A} = (vec{Q}, -vec{R})
end{align}

we can represent the constraint as:
begin{align}
vec{Ax} = vec{0}~~~textbf{ Constraint 1}
end{align}

The call prices at the input marks $hat k_j, j in [1, 2, ldots, M]$ can be represented by cubic spline interpolation
begin{align}
C(hat k_j) =a ci + b c{i+1} + frac{(a^3 a)u^2}6 pi + frac{(b^3-b) u^2}6 p{i+1}
~~~textbf{ Constraint 2}
end{align}

where

begin{align}
a = frac{k_{i+1} hat k_j}{u},~~~b = 1-a
end{align}

and $[k_i, k_{i+1}]$ here represents the segment that $hat k_j$ falls in.

$p_i$ are densities, so

begin{align}
p_i > 0 ~~~textbf{ Constraint 3}
end{align}

Integrating the density function we should get 1.0 (recall that density function are linearly interpolated)

begin{align}
u sum p_i = 1.0 ~~~textbf{ Constraint 4}
end{align}

Natural cubic spline, $p_1$ and $p_N$ are zero, so we could solve directly $c_1$ and $c_N$

begin{align}
c_1 = Se^{(r_d r_f)T} k_1, ~c_N = 0~textbf{ Constraint 5}
end{align}

Call prices are monotonically decreasing:

begin{align}
c{i+1} c{i} leq 0 ~text{for}~i in {1, 2, ldots, N-1}~textbf{ Constraint 6}
end{align}

Objective Function

Fill the rest of the DOF using objective function (soft constraints)

[Fengler 2009] tried minimizing the below to achieve smoothness on $p$:

begin{align}
int_{k_1}^{k_N} p(S)^2 dS = text{constant} times vec{p}^{top} vec{R} vec{p}
end{align}

Using $vec{x}$ as variable and define

begin{align}
vec{H}{(2N-2) times (2N-2)} =
begin{pmatrix}
vec{0} & vec{0}
vec{0} & vec{R}{(N-2) times (N-2)}
end{pmatrix}
end{align}

the problem becomes minimizing

begin{align}
vec{x}^top vec{H} vec{x}
end{align}

Problem Formulation

We can formulate our problem as
begin{align}
min~~~vec{x}^top vec{H} vec{x}
end{align}
subject to constraints 1 to 5.

All the constraints are linear function of $vec{x}$
Our objective function is quadratic and the matrix $vec{H}$ is positive semi-definite
Global solution exists, and (relatively) efficient to solve

Tips

To solve the quadratic programming problem, we can use the CVXOPT package:http://cvxopt.org/examples/tutorial/qp.html
https://courses.csail.mit.edu/6.867/wiki/images/a/a7/Qp-cvxopt.pdf

Write down the exact formulas using the same symbols used by CVXOPT QP problems documentation in the above docs, then translate them into code. This will make debugging easier.
To check whether solvers result makes sense, examine if the constraints are satisified, and if the call prices are smooth and match the input.
If test run takes too long, reduce the number of grid points in PDE pricer, or skip the calibration report and inspect the volatility surface first.
It might be easier to plot implied vol, call prices, PDF, and the marks to check the result.
use bisect.bisect_left to find the bucket $hat{k}$ belongs to (https://docs.python.org/3/library/bisect.html)

References

[Fengler 2009] Arbitrage-free smoothing of the implied volatility surface, Quantitative Finance, 2009

Implementation

Below are some building blocks for the project. You contribution should be in SmileAF class.
You can modify any other classes or methods. If you do so, please describe your modification in the project report.

Below are the smile interpolators and smile constructor. You need to implement SmileAF. Note that smileFromMarks takes a parameter smileInterpMethod. When it is AF, SmileAF is used.

Below is a calibration report that shows the calibration error of local volatility PDE pricer.

We test with no smile case first. In the calibration error report, we are showing the error in basis points 0.01% with respect to 1 notional.

In[4]:

S, r, q = 1.25805, 0.01, 0.0
iv = createTestImpliedVol(S, r, q, sc = 0.0, smileInterpMethod=CUBICSPLINE)
plotTestImpliedVolSurface(iv)
pdeCalibReport(S, r, q, iv)

Then test smile case with CubicSpline, with a mild smile (tuned by the coeffiicent sc)

In[5]:

iv = createTestImpliedVol(S, r, q, sc = 0.5, smileInterpMethod=CUBICSPLINE)
plotTestImpliedVolSurface(iv)
pdeCalibReport(S, r, q, iv)

Then test smile case with CubicSpline, with a the input smile (sc = 1.0). It can be seen that the short end low strike region has some smile arbitrage. The calibration errors become larger.

In[6]:

iv = createTestImpliedVol(S, r, q, sc = 1.0, smileInterpMethod=CUBICSPLINE)
plotTestImpliedVolSurface(iv)
pdeCalibReport(S, r, q, iv)

Your test cases with SmileAF

In[7]:

iv = createTestImpliedVol(S, r, q, sc = 1.0, smileInterpMethod=AF)
plotTestImpliedVolSurface(iv)
pdeCalibReport(S, r, q, iv)

ValueErrorTraceback (most recent call last)
in
-> 1 iv = createTestImpliedVol(S, r, q, sc = 1.0, smileInterpMethod=AF)
2 plotTestImpliedVolSurface(iv)
3 pdeCalibReport(S, r, q, iv)

in createTestImpliedVol(S, r, q, sc, smileInterpMethod)
19 bf10s = [0.0050, 0.0050, 0.0067, 0.0088, 0.0111, 0.0144, 0.0190, 0.0201, 0.0204, 0.0190, 0.0186, 0.0172, 0.0155, 0.0148]
20 rr10s = [-0.0111, -0.0187, -0.0248, -0.0315, -0.0439, -0.0518, -0.0627, -0.0652, -0.0662, -0.0646, -0.0636, -0.0627, -0.0627]
> 21 smiles = [smileFromMarks(pillars[i], S, r, q, atmvols[i], bf25s[i]*sc, rr25s[i]*sc, bf10s[i]*sc, rr10s[i]*sc, smileInterpMethod) for i in range(len(pillars))]
22 return ImpliedVol(pillars, smiles)
23

in (.0)
19 bf10s = [0.0050, 0.0050, 0.0067, 0.0088, 0.0111, 0.0144, 0.0190, 0.0201, 0.0204, 0.0190, 0.0186, 0.0172, 0.0155, 0.0148]
20 rr10s = [-0.0111, -0.0187, -0.0248, -0.0315, -0.0439, -0.0518, -0.0627, -0.0652, -0.0662, -0.0646, -0.0636, -0.0627, -0.0627]
> 21 smiles = [smileFromMarks(pillars[i], S, r, q, atmvols[i], bf25s[i]*sc, rr25s[i]*sc, bf10s[i]*sc, rr10s[i]*sc, smileInterpMethod) for i in range(len(pillars))]
22 return ImpliedVol(pillars, smiles)
23

in smileFromMarks(T, S, r, q, atmvol, bf25, rr25, bf10, rr10, smileInterpMethod)
100 return SmileCubicSpline(ks, [p10, p25, atmvol, c25, c10])
101 elif smileInterpMethod == AF:
> 102 return SmileAF(ks, [p10, p25, atmvol, c25, c10], T)

in __init__(self, strikes, vols, T)
38 f = lambda v: implyVol(self.ks[i], prc, v)
39 a, b = 1e-8, 10
> 40 vs[i (khmin-1) + 1] = optimize.brentq(f, a, b)
41 kks[i (khmin-1) + 1] = self.ks[i]
42 kks[0] = kmin

C:pathonlibsite-packagesscipyoptimizezeros.py in brentq(f, a, b, args, xtol, rtol, maxiter, full_output, disp)
774 if rtol < _rtol:775 raise ValueError(“rtol too small (%g < %g)” % (rtol, _rtol))–> 776 r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
777 return results_c(full_output, r)
778

ValueError: f(a) and f(b) must have different signs

In[]:

In[]: