The Black-Scholes equation is a partial differential equation that falls a bit out of the ordinary scheme. It describes what the fair price of a "European
call" stock option is. Without going into too much detail, a stock "option" is a contract one can buy from a bank that allows me, but not requires me, to buy a specific stock at a fixed price at a fixed future time in the future. The question one would then want to answer as a buyer of such an option is "How much do I think such a contract is worth?", or as the seller "How much do I need to charge for this contract?", both as a function of the time before the contract is up at time and as a function of the stock price . Fischer Black and Myron Scholes derived a partial differential equation for the fair price for such options under the assumption that stock prices exhibit random price fluctuations with a given level of "volatility" plus a background exponential price increase (which one can think of as the inflation rate that simply devalues all money over time). For their work, Black and Scholes received the Nobel Prize in Economic Sciences in 1997, making this the first tutorial program dealing with a problem for which someone has gotten a Nobel Prize [black1973pricing].
The equation reads as follows:
where
The way we should interpret this equation is that it is a time-dependent partial differential equation of one "space" variable as the price of the stock, and is the price of the option at time if the stock price at that time were .
Particularities of the equation system
There are a number of oddities in this equation that are worth discussing before moving on to its numerical solution. First, the "spatial" domain is unbounded, and thus can be unbounded in value. This is because there may be a practical upper bound for stock prices, but not a conceptual one. The boundary conditions as can then be interpreted as follows: What is the value of an option that allows me to buy a stock at price if the stock price (today or at time ) is ? One would expect that it is plus some adjustment for inflation, or, if we really truly consider huge values of , we can neglect and arrive at the statement that the boundary values at the infinite boundary should be of the form as stated above.
In practice, for us to use a finite element method to solve this, we are going to need to bound . Since this equation describes prices, and it doesn't make sense to talk about prices being negative, we will set the lower bound of to be 0. Then, for an upper bound, we will choose a very large number, one that is not very likely to ever get to. We will call this . So, .
Second, after truncating the domain, we need to ask what boundary values we should pose at this now finite boundary. To take care of this, we use "put-call" parity [stoll1969relationship]. A "pull option" is one in which we are allowed, but not required, to sell a stock at price to someone at a future time . This says
where is the value of the call option, and is the value of the put option. Since we expect as , this says
and we can use this as a reasonable boundary condition at our finite point .
The second complication of the Block-Scholes equation is that we are given a final condition, and not an initial condition. This is because we know what the option is worth at time : If the stock price at is , then we have no incentive to use our option of buying a price because we can buy that stock for cheaper on the open market. So for . On the other hand, if at time we have , then we can buy the stock at price via the option and immediately sell it again on the market for price , giving me a profit of . In other words, for . So, we only know values for at the end time but not the initial time – in fact, finding out what a fair price at the current time (conventionally taken to be ) is what solving these equations is all about.
This means that this is not an equation that is posed going forward in time, but in fact going backward in time. Thus it makes sense to solve this problem in reverse by making the change of variables where now denotes the time before the strike time .
With all of this, we finally end up with the following problem:
Conceptually, this is an advection-diffusion-reaction problem for the variable : There is both a second-order derivative diffusion term, a first-order derivative advection term, and a zeroth-order reaction term. We can expect this problem to be a little bit forgiving in practice because for realistic values of the coefficients, it is diffusive dominated. But, because of the advective terms in the problem, we will have to be careful with mesh refinement and time step choice. There is also the issue that the diffusion term is written in a non-conservative form and so integration by parts is not immediately obvious. This will be discussed in the next section.
Scheme for the numerical solution
We will solve this problem using an IMEX method. In particular, we first discretize in time with the theta method and will later pick different values of theta for the advective and diffusive terms. Let approximate :
Here, is the time step size. Given this time discretization, we can proceed to discretize space by multiplying with test functions and then integrating by parts. Because there are some interesting details in this due to the advective and non-advective terms in this equation, this process will be explained in detail.
So, we begin by multiplying by test functions, :
As usual, (1) becomes and (4) becomes , where , and where we have taken the liberty of denoting by not only the function but also the vector of nodal values after discretization.
The interesting parts come from (2) and (3).
For (2), we have:
There are two integrals here, that are more or less the same, with the differences being a slightly different coefficient in front of the integral, and a different time step for V. Therefore, we will outline this integral in the general case, and account for the differences at the end. So, consider the general integral, which we will solve using integration by parts:
So, after adding in the constants and exchanging for where applicable, we arrive at the following for (2):
But, because the matrix involves an advective term, we will choose there – in other words, we use an explicit Euler method to treat advection. Conversely, since the matrix involves the diffusive term, we will choose there – i.e., we treat diffusion using the second order Crank-Nicolson method.
So, we arrive at the following:
Now, to handle (3). For this, we will again proceed by considering the general case like above.
So, again after adding in the constants and exchanging for where applicable, we arrive at the following for (3):
Just as before, we will use for the matrix and for the matrix . So, we arrive at the following for (3):
Now, putting everything together, we obtain the following discrete form for the Black-Scholes Equation:
So, altogether we have:
As usual, we can write this with the unknown quantities on the left and the known ones on the right. This leads to the following linear system that would have to be solved in each time step:
Test Case
For this program, we will use the Method of Manufactured Solutions (MMS) to test that it is working correctly. This means that we will choose our solution to be a certain function similar to step-7. For our case, we will use:
This means that, using our PDE, we arrive at the following problem:
Where, . This set-up now has right hand sides for the equation itself and for the boundary conditions at that we did not have before, along with "final" conditions (or, with -time "initial conditions") that do not match the real situation. We will implement this in such a way in the code that it is easy to exchange – the introduction of the changes above is just meant to enable the use of a manufactured solution.
If the program is working correctly, then it should produce (**) as the solution. This does mean that we need to modify our variational form somewhat to account for the non-zero right hand side.
First, we define the following:
So, we arrive at the new equation:
We then solve this equation as outlined above.
The commented program
Include files
The program starts with the usual include files, all of which you should have seen before by now:
Then the usual placing of all content of this program into a namespace and the importation of the deal.II namespace into the one we will work in. We also define an identifier to allow for the MMS code to be run when MMS is defined. Otherwise, the program solves the original problem:
This section creates a class for the known solution when testing using the MMS. Here we are using for the solution. We need to include the solution equation and the gradient for the H1 seminorm calculation.
In the following classes and functions, we implement the right hand side and boundary values that define this problem and for which we need function objects. The right hand side is chosen as discussed at the end of the introduction.
The next piece is the declaration of the main class of this program. This is very similar to the step-26 tutorial, with some modifications. New matrices had to be added to calculate the A and B matrices, as well as the vector mentioned in the introduction. We also define the parameters used in the problem.
maximum_stock_price: The imposed upper bound on the spatial domain. This is the maximum allowed stock price.
maturity_time: The upper bound on the time domain. This is when the option expires.
asset_volatility: The volatility of the stock price.
interest_rate: The risk free interest rate.
strike_price: The agreed upon price that the buyer will have the option of purchasing the stocks at the expiration time.
Some slight differences between this program and step-26 are the creation of the a_matrix and the b_matrix, which is described in the introduction. We then also need to store the current time, the size of the time step, and the number of the current time step. Next, we will store the output into a DataOutStack variable because we will be layering the solution at each time on top of one another to create the solution manifold. Then, we have a variable that stores the current cycle and number of cycles that we will run when calculating the solution. The cycle is one full solution calculation given a mesh. We refine the mesh once in between each cycle to exhibit the convergence properties of our program. Finally, we store the convergence data into a convergence table.
As far as member functions are concerned, we have a function that calculates the convergence information for each cycle, called process_solution. This is just like what is done in step-7.
Now, we get to the implementation of the main class. We will set the values for the various parameters used in the problem. These were chosen because they are fairly normal values for these parameters. Although the stock price has no upper bound in reality (it is in fact infinite), we impose an upper bound that is twice the strike price. This is a somewhat arbitrary choice to be twice the strike price, but it is large enough to see the interesting parts of the solution.
The next function sets up the DoFHandler object, computes the constraints, and sets the linear algebra objects to their correct sizes. We also compute the mass matrix here by calling a function from the library. We will compute the other 3 matrices next, because these need to be computed 'by hand'.
Note, that the time step is initialized here because the maturity time was needed to compute the time step.
Below is the code to create the Laplace matrix with non-constant coefficients. This corresponds to the matrix D in the introduction. This non-constant coefficient is represented in the current_coefficient variable.
Now we will create the A matrix. Below is the code to create the matrix A as discussed in the introduction. The non constant coefficient is again represented in the current_coefficient variable.
for (constauto &cell : dof_handler.active_cell_iterators())
{
cell_matrix = 0.;
fe_values.reinit(cell);
for (constunsignedint q_index : fe_values.quadrature_point_indices())
Finally we will create the matrix B. Below is the code to create the matrix B as discussed in the introduction. The non constant coefficient is again represented in the current_coefficient variable.
for (constauto &cell : dof_handler.active_cell_iterators())
{
cell_matrix = 0.;
fe_values.reinit(cell);
for (constunsignedint q_index : fe_values.quadrature_point_indices())
The next function is the one that solves the actual linear system for a single time step. The only interesting thing here is that the matrices we have built are symmetric positive definite, so we can use the conjugate gradient method.
This is simply the function to stitch the solution pieces together. For this, we create a new layer at each time, and then add the solution vector for that timestep. The function then stitches this together with the old solutions using 'build_patches'.
It is somewhat unnecessary to have a function for the global refinement that we do. The reason for the function is to allow for the possibility of an adaptive refinement later.
This next part is building the convergence and error tables. By this, we need to set the settings for how to output the data that was calculated during BlackScholes::process_solution. First, we will create the headings and set up the cells properly. During this, we will also prescribe the precision of our results. Then we will write the calculated errors based on the , , and norms to the console and to the error LaTeX file.
Now we get to the main driver of the program. This is where we do all the work of looping through the time steps and calculating the solution vector each time. Here at the top, we set the initial refinement value and then create a mesh. Then we refine this mesh once. Next, we set up the data_out_stack object to store our solution. Finally, we start a for loop to loop through the cycles. This lets us recalculate a solution for each successive mesh refinement. At the beginning of each iteration, we need to reset the time and time step number. We introduce an if statement to accomplish this because we don't want to do this on the first iteration.
Next, we run the main loop which runs until we exceed the maturity time. We first compute the right hand side of the equation, which is described in the introduction. Recall that it contains the term . We put these terms into the variable system_rhs, with the help of a temporary vector:
vmult_result.reinit(dof_handler.n_dofs());
forcing_terms.reinit(dof_handler.n_dofs());
for (unsignedint timestep_number = 0; timestep_number < n_time_steps;
++timestep_number)
{
time += time_step;
if (timestep_number % 1000 == 0)
std::cout << "Time step " << timestep_number << " at t=" << time
The second piece is to compute the contributions of the source terms. This corresponds to the term . The following code calls VectorTools::create_right_hand_side to compute the vectors , where we set the time of the right hand side (source) function before we evaluate it. The result of this all ends up in the forcing_terms variable:
Next, we add the forcing terms to the ones that come from the time stepping, and also build the matrix that we have to invert in each time step. The final piece of these operations is to eliminate hanging node constrained degrees of freedom from the linear system:
There is one more operation we need to do before we can solve it: boundary values. To this end, we create a boundary value object, set the proper time to the one of the current time step, and evaluate it as we have done many times before. The result is used to also set the correct boundary values in the linear system:
Having made it this far, there is, again, nothing much to discuss for the main function of this program: it looks like all such functions since step-6.
What is more interesting is the output of the convergence tables. They are outputted into the console, as well into a LaTeX file. The convergence tables are shown above. Here, you can see that the solution has a convergence rate of with respect to the -norm, and the solution has a convergence rate of with respect to the -norm.