When constructing boundary conditions $\phi_1$ and $\phi_2$ in the case of option pricing it is often a good idea to use the fact, that the solution of the Black-Scholes equation satisfying the final condition $v(s,T)=c_1 +c_2 s$ is $$v(s,t)=c_1 e^{-r(T-t)}+c_2e^{-D(T-t)}s.$$ Hence, if the payoff function is a linear function starting from some value $s=s_1$, then for large values of $s$ the solution is practically equal to the special solution corresponding to this linear function. And if the payoff function is linear for $s<s_2$, than for small values of $s$ the solution is practically equal to the special solution corresponding to that linear function. This gives us a possibility to define boundary values so that they are not very different from the actual option prices at those boundaries and hence to reduce significantly the error caused by introducing $x_{min}$ and $x_{max}$ in option pricing equations. For example, when finding the value of a call option price by solving the untransformed equation, we have that for large values of $s$ the payoff is $p(s)=s-E$. Since for transformed equation $u(x,t)=v(e^x,t)$, a suitable boundary condition for $x=x_{max}$ is $$\phi_2(t, x_{max})= -E e^{-r(T-t)}+e^{-D(T-t)}e^{x_{max}}.$$
In the exercises 1 and 2 we consider an European option with $T=0.5$ and payoff function $$p(s)=\begin{cases} \frac{195-2\cdot s}{8},& s<95,\\ \frac{(s-100)^2}{40},& 95\leq s\leq 120,\\ s-110,& s>120. \end{cases}$$
Derive suitable boundary conditions (based on special solutions of BS equation) for the transformed Black-Scholes equation for the payoff function defined above. Define corresponding functions phi_1(t,xmin) and phi_2(t,xmax).
Solution
import numpy as np
def phi1(t,xmin):
c1=??
c2=??
return ??
def phi2(t,xmax):
c1=??
c2=??
return ??
If we have variable coefficient $\alpha(x,t)$ in the transformed Black-Scholes equation, then we can choose $m$ so that the stability constraint $b_{ik}\geq 0$ is satisfied for all $i,k$ if we know the maximal value of $\alpha_{max}=\alpha(x_i,t)$: $$b_{ik}=1-2\frac{T}{m\Delta x^2}\alpha(x_i,k\frac{T}{m})-r\frac{T}{m}\geq 1-2\frac{T}{m\Delta x^2}\alpha_{max}-r\frac{T}{m},$$ so by choosing m that satisfies the condition $$1-2\frac{T}{m\Delta x^2}\alpha_{max}-r\frac{T}{m}\geq 0$$ we make sure that the stability constraint is satisfied. If we want that our program works with all reasonable functions $\alpha$, we have to write a piece of code that estimates the value $\alpha_{max}$ automatically. A possible way to do it is to use the following steps
After estimating $\alpha_{max}$ we can find $m$ form the stability constraint and compute the approximate values of the solution of the partial differential equation by the explicit finite difference method (see the handout of the previous lab).
Define a function f_max_2d(f,x,T,t0=0,m0=100,c=10) which estimates the maximal value of a function of two variables when the first argument takes the values in the given vector $x$ and the second argument is between $t_0$ and $T$, by computing the values at $m_0+1$ points in the interval $[t_0,T]$ for each value of of the first argument from the vector $x$. Check the correctness of your function when $f(x,t)=\frac{10}{(x-1)^2+t+0.1}$, x is a uniformly spaced vector in the interval $[0,\,3]$ with 31 values and $T=1$ (other parameters use default values).
def f_max(f,x,T,t0=0,c=10,m0=100):
#x is a vector of points x_i
#c is an estimate to max value of derivative of f w.r.t t
#f is a function of x and t
return ???
Check the correctness of the function:
x=np.linspace(0,3,31)
T=1
def f(x,t):
return ??
print(f_max(f,x,T))
Implement the explicit finite difference method described in the previous lab so that it works with non-constant $\alpha$ and $\beta$. Test your solver by using the boundary conditions derived in the previous exercise in the case $n=20$, $\rho=1.5$, $r=0.01$, $D=0.05$, $T=0.5$, $S0=97$ and $$\sigma(s,t)=0.5+0.2\cdot e^{-t}\,\arctan(0.1\,s-10).$$ The answer depends on how a suitable value of $m$ is determined from the stability condition, but should be close to 10.67.
Solution
Use the function f_max defined in the previous exercise for estimating the maximal value of alpha.
def explicit_solver2(n,rho,r,D,S0,T,sigma,p,phi1,phi2):
"""sigma is assumed to be a function of s and t
phi1,phi2 are functions of t,xmin and t,xmax
p is a function of stock price only
Requires that the function f_max which estimates the maximum value of a
function of two variables is defined.
"""
xmax=np.log(S0*rho)
xmin=np.log(S0/rho)
delta_x=(xmax-xmin)/n
x=np.linspace(xmin,xmax,n+1)
#find m from the stability condition
#for this, we first have to estimate the maximum of the coefficient alpha
#for transformed BS equation, alpha is defined as follows
def alpha(x,t):
return ??
#estimate max value of alpha
alpha_max=??
#Estimate m from the stability condition
m=??
m=np.int64(np.ceil(m)) #has to be an integer
delta_t=??
#continue similarly to lab9. Note that a, b, c are now vectors and should be computed
#inside for cycle for each time moment separately
return U[n//2,0] #option price at t=0 when S(0)=S0
Check the solver:
#define parameters
#define payoff function and sigma
def p(s):
return (195-2*s)/8*(s < 95)+??*(s >= 95)*(s <= 120) + ??
def sigma(s,t):
return ??
approx_sol=explicit_solver2(??)
print(approx_sol) #shoulb be approximately 10.67
Consider the European option with the pay-off function $$p(s)=\begin{cases} \frac{s}{4},& 0\leq s<40,\\ \frac{1}{100} \, {\left(s - 40\right)}^{3} - \frac{1}{6} \, {\left(s - 40\right)}^{2} + \frac{s}{4},& 40\leq s\leq 55,\\ 2\,(s-50)& s>55. \end{cases}$$ Assume $r=0.01$, $D=0.05$, $S(0)=50$, $T=1$ and that the volatility function is given by $$\sigma(s,t)=0.5+\frac{0.2\cdot e^{-t}}{1+0.01\cdot(s-40)^2}.$$ For $\rho=2.5$ and both in the case of the boundary conditions given by appropriate special functions and in the case of using the simple (constant) boundary conditions, determine the limiting value (as $n\rightarrow\infty$) of the approximate option price at time $t=0$ with the accuracy 0.005 by computing the approximate values by the explicit finite difference method for $n=10,20,40, ...$ and by estimating the accuracy of the last result by $$\frac{1}{3}|\,\text{lastResult}-\text{previousResult}\,|.$$ Repeat the computations in the case $\rho=8$ and accuracy 0.001. Which method gave more accurate results for $\rho=2.5$?
Solution
#define parameters and functions needed for the solver (boundary conditions,)
First, consider the case of constant boundary conditions
#compute the limiting value for the option price
Now the same with the special boundary conditions:
??
Computations for $rho=8$, first for special boundary conditions.
rho=8
limitingError=0.001
Repeat the same computations with the constant boundary conditions:
rho=8
limitingError=0.001