If we consider an European option with exercise time $T$ and payoff function $p$ and assume the validity of Black-Scholes market model, then the option price at time $t$ is given by $v(S(t),t)=u(\ln(S(t)),t)$, where $u$ is the solution of the problem \begin{equation*} \frac{\partial u}{\partial t}(x,t)+\alpha(x,t)\frac{\partial^2 u}{\partial x^2}(x,t)+ \beta(x,t)\frac{\partial u}{\partial x}(x,t)-r\, u(x,t)=0,\ x\in \mathbb{R}, 0\leq t<T \end{equation*} satisfying the final condition $$u(x,T)=p(e^x),\ x\in \mathbb{R}.$$ Here \begin{align*} \alpha(x,t)&=\frac{\sigma^2(e^x,t)}{2},\\ \beta(x,t)&=r-D-\frac{\sigma^2(e^x,t)}{2}. \end{align*} For solving the equation for $u$ numerically, we introduce two boundaries $x_{min}$ and $x_{max}$ and specify boundary conditions $u(x_{min},t)=\phi_1(t,x_{min}),\ u(x_{max},t)=\phi_2(t,x_{max})$ at those points. Next, we introduce the points $x_i=x_{min}+i\Delta x,\ i=0,\ldots,n$ and $t_k=k\Delta t,\ k=0,\ldots,m$ and denote by $U_{ik}$ approximate values of $u(x_i,t_k)$. Here $\Delta x=\frac{x_{max}-x_{min}}{n}$ and $\Delta t=\frac{T}{m}$. In the case of the explicit finite difference method we compute the values $U_{ik}$ as follows: \begin{align*} U_{im}&=p(e^{x_i}),\ i=0,\ldots,n\\ U_{0,k-1}&=\phi_1(t_{k-1}),\ U_{n,k-1}= \phi_2(t_{k-1}),\ k=m,m-1,\ldots,1,\\ U_{i,k-1}&=a_{ik}U_{i-1,k}+b_{ik}U_{ik}+c_{ik}U_{i+1,k},\ i=1,\ldots,n-1,\ k=m,m-1,\ldots,1, \end{align*} where \begin{align*} a_{ik}&=\frac{\Delta t}{\Delta x^2}\left(\alpha(x_i,t_k)-\frac{\beta(x_i,t_k)}{2}\Delta x\right),\\ b_{ik}&=1-2\frac{\Delta t}{\Delta x^2}\alpha(x_i,t_k)-r\Delta t,\\ c_{ik}&=\frac{\Delta t}{\Delta x^2}\left(\alpha(x_i,t_k)+\frac{\beta(x_i,t_k)}{2}\Delta x\right). \end{align*} If $\sigma$ satisfies $0<c_1 \leq \sigma(s,t)\leq c_2\ \forall s>0,\ \forall t\in [0,T)$ for some constants $c_1$ and $c_2$, then the method is stable for large enough $m,n$ if $b_{ik}\geq 0$. This means that we can choose one of the paramters $m,n$ freely, but the other one should be chosen so that the constraint is satisfied. Usually it is better to choose $n$ and then choose $m$ so that the condition $b_{ik}\geq 0$ is satisfied.
If $\sigma$ is a constant, then the coefficients $a=a_{ik},b=b_{ik}$ and $c=c_{ik}$ are also constants (do not depend on $i$ and $k$) and the numerical scheme simplifies to $$U_{i,k-1}=a\,U_{i-1,k}+b\,U_{ik}+c\,U_{i+1,k},\ i=1,\ldots,n-1,\ k=m,m-1,\ldots,1.$$ The stability condition is in this case $b\geq 0$.
If we write down this condition in terms of $m$ and $n$ we get $$1-2\frac{T}{m}\frac{n^2}{(x_{max}-x_{min})^2}\alpha-r\frac{T}{m}\geq 0,$$ where $\alpha=\frac{\sigma^2}{2}$.
Solve the stability condition for $m$, that is, find an equivalent condition of the form $m\geq \ldots$.
Solution: ???
Write a function that for given values of $n$, $\rho>1$, $r$, $D$, $S0$, $T$, $\sigma$ and for given functions $p$, $\phi_1$ and $\phi_2$ takes $m$ to be equal to the smallest integer satisfying the stability constraint (hint: use commands np.int64(np.ceil(...))) and returns the values $U_{i0},\ i=0,\ldots,n$ of the approximate solution (option prices) obtained by solving the tarnsformed BS equation with the explicit finite difference method and the corresponding stock prices $S_i=e^{x_i}$ in the case $x_{min}=\ln \frac{S0}{\rho},\ x_{max}=\ln(\rho\, S0)$. Test the correctness of your code by comparing the results to the exact values obtained by Black-Scholes formula in the case $r=0.03$, $\sigma=0.5$, $D=0.05$, $T=0.5$, $E=97$, $S0=100$, $p(s)=\max(s-E,0)$, $\phi_1(t)=p(e^{x_{min}}),\ \phi_2(t)=p(e^{x_{max}})$.
Solution
import numpy as np
def explicit_solver(n,rho,r,D,S0,T,sigma,p,phi1,phi2):
"""assume that phi1 and phi2 are functions of 2 arguments (xmin,t or xmax,t)"""
xmax=
xmin=
delta_x=
#find m from the stability condition
m=
#change it to an integer, which is not smaller than the computed value
m=
delta_t=
#define values of x_i
???
for k in range(m,0,-1): #backward iteration, k=m,m-1,...
#boundary conditions
#all other values
return [U[:,0],np.exp(x)]
Test the solution. For this, define the pay-off function and the functions $\phi_1,\ \phi_2$
#define parameters
#define payoff function
def p_call(s):
return
#define boundary functions
def phi1(xmin,t):
return
def phi2(xmax,t):
return
n=6;rho=4
solution=explicit_solver(n,rho,r,D,S0,T,sigma,p_call,phi1,phi2)
print("option prices:",solution[0])
print("stock prices:",solution[1])
Note that the relation to the matrix U computed inside the function explicit_solver and the option prices are as follows:
if at the time moment $t=t_k$ the logarithm of the current stock price $\ln(S(t))$ is $x_i$, then the price is approximately $U_{ik}$. Since we return only to values corresponding to $t=t_0=0$, the first vector of the result gives the prices for different possible values of $S(0)$: if $\ln(S(0))=x_i$, then the option price is the $i$-th value in the first vector.
So $x_{min},\ x_{max}$ and $n$ determine, for which stock prices we have approximate values in the matrix $U$ (and thus also in the first component of the result of the function). Now, if we know in advance that we are mostly interested in the case $S(0)=S_0$, where $S_0$ is a given number, it is good to make sure that $\ln S_0$ is one of the grid points in $x$ direction.
This is why we have defined $x_{min}$ and $x_{max}$ by the formulas given in the previous exercise handout. Namely, since $\frac{x_{min}+x_{max}}{2}=\ln S_0$ (show it!), we have the price corresponding to $S(0)=S_0$ in the matrix $U$ whenever $\frac{x_{min}+x_{max}}{2}$ is a grid point, and that is true for all even values of $n$. Moreover, the price is exactly in the middle of the $n+1$ values our function returns, which is clearly seen also by looking at the stock prices in the output: the value $S_0=100$ is exactly the middle element of the vector of stock prices.
Compare the approximate price of the option to the exact value computed by the BS formula:
n=6;rho=4
The difference is quite small even for $n=6$. Check the result for a larger value of $n$:
n=100
Let $r=0.02$, $\sigma=0.6$, $D=0.03$, $T=0.5$, $E=99$, $S0=100$, $p(s)=\max(s-E,0)$. If we use the explicit method of previous exercise, then even if we let $n$ go to infinity there is going to be a finite error between the exact option price at $t=0, S(0)=S0$ and the corresponding approximate value. This error is caused by introducing artificial boundaries $x_{min}$ and $x_{max}$ and the boundary conditions specified at those boundaries. Use the boundary conditions $\phi_1(t)=p(e^{x_{min}}),\ \phi_2(t)=p(e^{x_{max}})$ and determine the value of the resulting error for $\rho=1.5, 2, 2.5$. In order to see the resulting error you should do several computations with fixed $\rho$ and increasing values of $n$ (assuming $m$ is determined from the stability condition, $n$ should be increased by multiplying it by 2 each time). Use the knowledge that for large enough $n$ the part of the error depending on the choice of $n$ behaves approximately like $\frac{const.}{n^2}$ (so the difference of the last two computations divided by 3 is an estimate of this part of the error for the last computation) for determining how far your last computation is from the limiting value. For each value of $\rho$, start computations from $n=10$ and stop the computations when the absolute error depending on $n$ is at least 10 times smaller than the actual absolute error of the computed answer.
Solution
#define parameters, compute the exact price
First value of $\rho$:
rho=1.5
n=10
The second value of $\rho$:
rho=2
n=10
So the limiting error is approximately 0.034
The last value of $\rho$:
rho=2.5
n=10
So we see that the error caused by the location of $x_{min}$ and $x_{max}$ is approaching 0 quickly when $\rho$ increases. For most practical problems $\rho=3$ is good enough but of cause we should do some computations to check if this is indeed the case.
It is important to understand that a single computation by a finite difference method can be used for computing prices of the same option at different time moments for any reasonable value of the current stock price at the time moments of interest. Of cause our numerical method does not give us a formula for the function $v$, but only a table of values corresponding to certain stock prices and time moments. So to find a value of the function for a given value of $s$ and for given time moment, we should determine the closest values of of stock prices and time moments for which we have results and to use them to compute the value we need. One commonly used approach is bilinear interpolation, which works as follows.
Suppose that we know the values of a function $u(x,t)$ at the corners of a rectangle $[x_1,x_2]\times [t_1,t_2]$, denote the known values by $u(x_1,t_1)=U_{11}$, $u(x_1,t_2)=U_{12}$, $u(x_2,t_1)=U_{21}$, $u(x_2,t_2)=U_{22}$. Denote also $\Delta x=x_2-x_1,\Delta t=t_2-t_1$. Then the the value of $u$ at the point $(x,t)$, where $x_1\leq x<x_1$, $t_1\leq t\leq t_2$ can be computed approximately by the bilinear interpolation formula $$u(x,t)\approx \frac{U_{11}(x_2-x)(t_2-t)+U_{12}(x_2-x)(t-t_1)+U_{21}(x-x_1)(t_2-t)+U_{22}(x-x_1)(t-t_1)}{\Delta x\,\Delta t}.$$ It is known that if the function $u$ has continuous second derivatives in $x$ and $t$ variables, then the approximation error is bounded by $const. \cdot (\Delta x^2+\Delta t^2)$.
Modify the solver of Exercise 2 so that in addition of $S0$ and $\rho$ it takes the value of $x_{min}$ as an argument of the function, uses $x_{max}=\ln(\rho\,S0)$ for the location of the upper boundary and returns the matrix $U$ as the answer. Use the explicit finite difference method with $r=0.01$, $\sigma=0.6$, $D=0$, $T=0.7$, $S0=38$, $p(s)=\frac{30}{1+0.3(s-40)^{2}}$, $x_{min}=3,\ \rho=2$, $\phi_1(t,x_{min})=p(e^{x_{min}}),\ \phi_2(t,x_{max})=p(e^{x_{max}})$, $n=50$ to compute the matrix $U$. The matrix should be computed only once when solving the homework problem. Use this matrix to find approximate option prices:
I ask you to solve the problem by yourself without discussing it with other students(you can ask for help from me)!