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}
or
\begin{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 measure
\begin{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) ds
\end{align}The 1st non-arbitrage condition translates to
\begin{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{K\int_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 -\frac12\sigma_{ATM}^2T – 5 \sigma_{ATM} \sqrt{T}}, k_N = S e^{(r_d – r_f)T -\frac12\sigma_{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 first derivative 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 problem’s documentation in the above docs, then translate them into code. This will make debugging easier.
To check whether solver’s 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)
—————————————————————————
ValueError Traceback (most recent call last)
—-> 1 iv = createTestImpliedVol(S, r, q, sc = 1.0, smileInterpMethod=’AF’)
2 plotTestImpliedVolSurface(iv)
3 pdeCalibReport(S, r, q, iv)
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
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
100 return SmileCubicSpline(ks, [p10, p25, atmvol, c25, c10])
101 elif smileInterpMethod == “AF”:
–> 102 return SmileAF(ks, [p10, p25, atmvol, c25, c10], 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:\pathon\lib\site-packages\scipy\optimize\zeros.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 [ ]: