Extreme Vasicek is not Extreme Enough

Recently, in Extreme Vasicek, I showed some results, when the (absolute) volatility of a Vasicek model was by a factor 50 or so higher than it occurs in reality.

I confess: I was so satisfied with the quality of our results that I did not ask the correct question: Is this high volatility really a difficult case?

High volatility adds diffusion to the differential equation which should be then easy to be solved as long as you use implicit or semi-implicit schemes. Standard finite element or fintie difference schemes (with central differences used as approximation for the first order derivatives) might get into troubles, when convection gets dominant. In the Vasicek exmaple, this means that volatility should be low, reversion speed should be high for the test case.

 Therefore, we now try
(1)    a = 0.05,    b = 0.5,    sigma = 0.005

and

(2)    a = 0.05,    b = 0.5,    sigma = 0.0005

and calculate the prices V of the 10 years zero coupon bond for an initial short rate r = 0.

The analytic solution gives
(1)  V= 0.4489
and
(2) V = 0.4487.

We use Dirichlet conditions V=1 at r = -0.4 and at r = 0.4 and compare to numerical schemes
(a) fully implicit time-stepping, central differences for first derivatives
(b) fully implicit time-stepping, upwinding for first derivatives.

Here are the results.

Case (1)
Method (a):   V = 0.4489
Method (b):   V = 0.4492

Case (2)
Method (a):   V = -0.0509
Method (b):   V = 0.4491.

What happened here?
The bond prices for different initial spot rates calculated by Method (a) exhibit the following behavior in case(2)

whereas the upwinding method (b) shows an almost boring plot

Why could we observe an obviously wrong result in case (1)? Well, we could have. It was just (good? bad?) luck that the oscillations did not spread to the center yet.

The upwinding method shows the same smooth behavior as before.

My conclusion when convection plays a role: Use stable methods to get reliable results.
By the way, you can find such methods in the Workout

Hows, Whys and Wherefores of FEM in Quantitative Finance - V

Elements II

In the last blog post we started discussing general properties of elements and took a look at a simple one dimensional element. Today we will analyse elements for a 2D setup for which we introduce the triangular element families. In both cases we focus on linear and quadratic shape functions, which we consider sufficient for financial modeling.

The two-dimensional linear triangular element is also known as simplex element

and is represented by a linear polynomial in x and y (the superscript (e) denotes an element):

The nodal values can again be used to determine the coefficients by solving a linear system (a little homework)

yielding

and an unknown quantity can again be expressed using the nodal values and the shape functions.

Omitting the superscript (e), the coefficients a, b, and c int the shape functions can be easily calculated using the node co-ordinates

Note that these shape functions satisfy the properties discussed for the one dimensional case: Each shape function has a value of one at its own node and a value of zero at the other two. Summing up all three functions always amounts to one. The shape functions vary linearly along the edges between its node and the other two nodes, and the shape function is zero along the edge opposite its node. Also note that derivatives with respect to the spatial coordinates are constant within an element.

In our book (see link below) we also present the rectangular element family and discuss numerical integration techniques used for higher order and isoparametric elements.

This series of blog posts summarizes a chapter of the book A Workout in Computational Finance

The sixth post of the series will be published on Thursday, July 4 and will show how the 
assembling of a linear system of equations using the discrete finite elements works out. We show this procedure in a step-by-step manner by applying the finite element method to a two dimensional static diffusion-reaction equation.

Extreme Vasicek Examples the World Was (not) Waiting for

My recent blog
Setting boundary conditions that you don't know

was followed by activated discussions on Wilmott's book forum.
Alan suggested to do similar calculations for more difficult examples than the one I used there, and this is exactly what I did.

Again, we work on the Vasicek model

but now sigma = 0.5 (that is 50 percent absolute volatility for the short rate), which is completely unrealistic for interest rates in the main currencies. The reversion speed should be b=0.5, and the parameter a either a=0 or a=1 (the second case means a mean reversion level of 200 percent). The bond values V for the 10 year zero coupon bonds with face amounts 1 can be calculated analytically (e.g., following Shreve). Assuming the random walks (in both cases) to start at r = 0, they are 

sigma = 0.5, b=0.5, a=0:   V = 33.56

and

sigma= 0.5, b=0.5, a=1:    V = 3.68 E-06.

Of course, as the volatility is much higher than in the past blog, the calculation domain has to be chosen larger. Here are the results for the first case:

For the artificial boundary conditions V = 1, set at r = - 4 and at r = 4, we obtain (with a mesh size of 0.002) a bond value of 33.35. Again, we used upwinding for stabilization, which was quite successful, as indicated by the smooth curve.

The discrepancy between the analytic value and the numerical value stems from the mesh size (20 basis points), not from the boundary conditions.

In the second case (mean reverting to 200 percent), we obtain for two different Dirichlet boundary conditions to be set at r = -2 and at r = 6 

Dirichlet boundary conditions V=0

Dirichlet boundary conditions V=1

In both cases, the numerical value for the bond was 3.48E-06 for a mesh size of 0.01 (1 percent).
The difference of the values for the two Dirichlet conditions can be seen (as a matter of fact, can not be seen) here:

No visible difference between the two Dirichlet boundary values outside the boundary layers

My conclusions are:
1) Upwinding works excellently, also for weird cases
2) Boundary conditions do not influence the instrument value as long as set reasonably far away from the reversion level.
3) Mesh size does influence the value due to numerical errors.

Why Is It That Socks Always Get Lost in The Laundry?

From Ask Ariely .. I really like his answer - in short: you have many socks, socks are in pairs and if you miss one and find it later elsewhere, but forgot that you missed it once, you miss its partner ... You count as lost each sock in a pair - even if non of them is really lost (in the laundry).

How to Do Lab and Factory Work Simultaneosly

Motivated by Seth Godin's post:

At the lab we strive for finding a breakthrough, new ways to solve new problems and do new things. In labs we accept to get insight from failure. A lab does not think too much about exploiting, but creating.

Setting Boundary Conditions that You Don’t Know

Assume that you are somehow involved in the aerodynamic imbalance estimates for wind turbines which begin to vibrate if not properly balanced. For doing a wind simulation of the rotors (comparable to a wind tunnel simulation for a car), you must restrict your calculation domain to a region in the neighborhood of the wind turbine. The air pressure and the wind speed at the boundary of your calculation domain are then set more or less arbitrary, and dangerous cases have to be found out systematically or by educated guesses. (The dangerous cases have to be un-risked. Sorry for that one.) By the way: Read Chapter 5 of the Workout to learn about upwinding, flowlines and appropriate discretisations of convection-dominated flows (and also why this is relevant for your finance application).

We study the influence of boundary conditions for a one-factor Hull-White model, here, for the sake of simplicity, a Vasicek model. The short rate r satisfies the stochastic differential equation

with the usual notation. The parameter b, the so-called reversion speed, pulls back the short rate to the equilibrium. This means that short rates far away from the equilibrium are very unlikely and that possibly wrong boundary conditions should not play a major role as long as these artificial boundaries are far away from the equilibrium. Can we confirm this argument?

We valuate a zero coupon bond that matures in 10 years from today. The parameters of the Vasicek model are chosen as b = 0.1, sigma = 0.01, and a = 0 meaning that negative interest rates can and will occur. We set artificial (Dirichlet) boundary conditions for the value of the bond at r = - 0.2 and at r = 0.2  and assume these boundary values to be 0, 2, 4 or 6 times the notional value. To solve the Vasicek differential equation, we use implicit time-stepping and appropriate upwinding to treat the convection correctly. Here are the present values of the bond  (for the different boundary values) as function of today’s short rate.

We observe that the boundary layers (the region influenced by the boundary conditions) are very thin. Of course, the calculation domain has to be chosen carefully: It should not be too small for accuracy reasons nor too big for performance reasons.

Read more about boundary conditions in Chapter 6 of the Workout

Hows, Whys and Wherefores of FEM in Quantitative Finance - IV

Elements I

In today's blog post we start giving an overview over different element and the shape functions N_i(x) used for generating the trial functions h(x). Independent of the dimension we first list a number of general properties valid for all shape functions we will discuss:

(1) A shape function is always one at its designated node (same nodal index) and zero on the other nodes of the element

(2) at any point in the element (including the boundaries) the sum over all shape functions is equal to one

Many problems in computational finance can be described by one-dimensional models, e.g., instruments under one factor short rate models or all options with only one underlying. The one-dimensional region is then a line and the division into elements or subregions is very straightforward. The concept of interpolation of a certain order between nodal values we show in the following for the one dimensional case can easily be applied across orders and dimensions.

1D linear Element:

A line element with two end nodes X_i and X_j (X_i < X_j) has the length L=X_j - X_i and corresponding nodal values Φ_i and Φ_j. Assuming a linear interpolation such that the unknown function Φ(x) varies linearly between the nodes within the actual element we obtain: 

 The coefficients α₁ and α₁ can be determined from the nodal values by solving:

Such equations, where the nodal values are multiplied by affine linear functions of x, are typical for the finite element formalism. The linear functions appearing in the equations above are called shape functions or interpolation functions. In the literature, the variable N with a subscript that indicates the node of the element is frequently used to denote shape functions. 

Often vector notation is used to display the equations.

The figure below displays the shape functions N_i and N_j:

It is easy to see, that the general properties listed at the beginning of the post are fulfilled by our shape functions.

This series of blog posts summarizes a chapter of the book A Workout in Computational Finance

The fifth post of the series will be published on Thursday, June 27 and will show how to apply the concept we derived in today's post can be extended to higher orders and higher dimensions.

UnRisk 7 Rolls Out

17-Jun-13 - we have released UnRisk PRICING ENGINE and UnRisk-Q version 7, introduced as UnRisk 7. This release is free for all UnRisk Premium Service Customers and will be shipped to all new customers immediately. UnRisk has been introduced 2001. Now UnRisk 7 is the 20th release.

Hows, Whys and Wherefores of FEM in Quantitative Finance - III

Grid Generation

In contrast to physics and engineering, where the optimal element type and approximation order is often determined by the nature of the problem, the situation is not equally clear when applying the finite element method to financial problems. Two key questions need to be asked in advance of tackling a problem:

(1) Is it necessary to dynamically adapt the mesh during the simulation process (for example by locally adding additional elements to refine the mesh)?

(2) Are derivatives of the (approximate) solution required with high order?

To answer the first question we look at the problems's spatial dimension. A simple Hull & White Model leading to a 1+1 (spatial and time dimension) dimensional partial differential equation (PDE) can be solved fast without the need of applying refinement techniques. This picture changes when for example a corss currency swap needs to be valuated. Using two one factor short rate models (for example two Hull & White models) and a one factor model (Garman and Kohlhagen) for the FX rate one will end with a 3+1 dimensional PDE. In this case a clever mesh refinement strategy can boost your calculation.

The answer to the second question indicates the required order of approximation: If the calculation of greeks or similar quantities connected with derivatives in higher order is of relevance, higher order approximations should be used.

Although for technical applications, the mesh generation algorithms employed, generally involve complex iterative smoothing techniques that attempt to align elements with boundaries or physical domains, in the field of financial engineering structured rectangular or hexahedral meshes can be easily generated using tensor products (cartesian grids) leading to structured meshes. Strictly speaking, a structured mesh, like the one in the figure below, can be identified by all interior nodes of the mesh having an equal number of adjacent elements. These meshes may also serve as starting points for structured triangular or tetrahedral meshes.

Structured Mesh

Unstructured meshes, on the other hand, drop the node valence requirement and allow any number of elements to meet at a single node. Triangular and tetrahedral meshes are by far the most common forms used here.

Unstructured mesh. Additionally local grid refinement is applied. 

Most mesh generation techniques for unstructured meshes currently in use fall into one of three categories, namely the Quad/Octree technique, Delauny triangulation and the advancing front algorithm.

The Octree technique recursively subdivides cubes or quads containing the computational domain until the desired resolution is achieved. Today methods based on Delauny criterion are by far the most popular methods among the triangular and tetrahedral meshing techniques. Put simply, the Delaunay criterion states that any node must not be contained inside the circumsphere/circumcircle of any tetrahedra/- triangle within the mesh. Another popular family of triangular and tetrahedral mesh generation algorithms is the advancing (or moving) front method. Tetrahedra are progressively built from the triangulated surface of the computational domain inward. An active front is maintained where new tetrahedra/triangles are formed.

To obtain unstructured quadrilateral or hexahedral meshes, unstructured triangular or tetrahedral meshes often serve as starting points. 

This series of blog posts summarizes a chapter of the book A Workout in Computational Finance

The fourth post of the series will be published on Thursday, June 20 and will give an overview over different element types.

Skateboarding and Computational Finance

My son is an avid skateboarder, and this is the reason why I know two or three of skateboarding vocabulary. I assume that the general public would name an obstacle like the following one as halfpipe but the correct word is miniramp because it (the miniramp) has no vertical boudaries.

The Miniramp Option was plotted as the value of a chooser option: At expiry of the chooser option, the holder of the chooser option has the right to decide if the chooser option should be exchanged into a call or into a put option on the same underlying with the same strike price at a future expiry date of the respective vanilla options.

Chapter 6 of the Workout is aimed at formulating terminal conditions, interface conditions (as the above menstioned chooser interface condition) and boundary conditions, which are of particular importance when the domain of calculation is, in  principle, unbounded. I will discuss this boundary conditions aspect in a future post.

Advanced skaters enjoy gaps and jumps. Therefore I introduced two discrete dividend payments into the miniramp option. Can you explain why the dividend steps are non-symmetric?

If you cannot explain it: You can preorder the Workout.

What Is It That Makes a Modern Product Manager?

I took a few days off and just returned from the Drome region in South France. I found a real quiet place, in a tiny village with 56 inhabitants, close to Crest with a stunning view into a green valley.

Hows, Whys and Wherefores of FEM in Quantitative Finance - II

Weighted residual methods
After the short introduction into the Finite Element Method in the previous post, today we want to examine the different weighted residual methods.
If an approximate solution (trial function) is substituted into a differential displaymath a residual remains, since the approximate solution does not satisfy the equation. Consider, as an example, the one dimensional second order partial differential equation

with appropriate initial/terminal/boundary conditions. Substituting y(x) by the approximative solution h(x),

yields a non-zero residual R(x) since h(x) does not satisfy the equation exactly. The weighted residual method demands that the integral over the domain over the residual times a weighting function is equal zero.

The number of weighting functions W_i(x) is determined by the number of unknown coefficients in the trial function. Different choices for the weighting functions --- which are sometimes also called testing functions --- are available. We list some of the more popular choices here:

(1) Collocation method:Dirac-delta functions are used as weighting functions,
W_i(x)=δ(x-X_i). The residual will thus vanish at the points X_i if the weighted residual equation is fulfilled.

(2) Subdomain method: the idea is to force the weighted residual to be zero not just at fixed points in the domain, but over various subsections of the domain. To accomplish this, the weight functions are set to unity (W(x)=1), and the integral over the entire domain is broken into a number of subdomains sufficient to evaluate all unknown parameters.

(3) Galerkin method: if we consider a trial function of the form

the Galerkin method uses the same functions N_i(x) for W_i(x) that are used for the approximate solutions. Under the assumption that we are looking for a function y in a function space U also the testing functions W are element of the function space U (for example polynomials of order one). This method is the most popular one and is used for solving problems with first derivative terms. We will focus on the Galerkin method exclusively throughout this blog post series.

(4) Least squares method: the residual is used as a weighting function and a new error term is defined which is then minimised with respect to the unknown coefficients of the trial functions.

This series of blog posts summarizes a chapter of the book A Workout in Computational Finance

The third post of the series will be published on Thursday, June 13 and will give an overview over different types of grids and grid generating algorithms.

Hows, Whys and Wherefores of FEM in Quantitative Finance - I

An Introduction to the Finite Element Method
This is the first part of a concise series of blog posts dealing with the Finite Element Method in quantitative finance. Today's blog is essentially a short introduction to the method and an outline of the basic steps necessary for application of the method to practical problems. Each step will be discussed in a blog entry.

The finite element Finite Elements Method is a numerical method for solving partial differential equations (PDEs), and has become particularly popular in engineering and physics. More recently also the quantitative finance community has gained interest to apply this method to the various PDEs arising in the modelling of financial instruments.

The basic strategy of the finite element method is to divide the region of interest into smaller parts (finite elements), and to approximate the solution in each element by a simple function. The method can thus be seen as a piecewise approximation. Polynomial-type interpolation functions are most widely used for the element solutions. Choosing polynomials over other kind of interpolation functions is convenient because of two reasons:

(1) It is relatively easy to formulate and compute the finite element equations.
(2) It is possible to obtain higher accuracy by increasing the order of the polynomial.

The practical application of the finite element method can be broken down into the following series of basic steps:

(1) Discretiztion of the region of interest: The most important concepts used for grid generation, such as refinement strategies and adaptivity, will be presented during this blog series. This step also includes numbering the nodes and specifying their corresponding coordinate values.
(2) The trial function h(x) is specified (for instance, the order of a polynomial approximation) and for each element, the partial differential equation is written in terms of the unknown nodal values.
(3) Using the Galerkin method (will be discussed in the next part of the blog series) for each element a system of equations is developed. These equations are then assembled into the global system of equations that covers the whole domain.
(4) Add the boundary/terminal/interface conditions
(5) Solve the system of equations for different methods to solve linear systems of equations

The second post of the series will be published on Friday, June 7 and will give an overview over different weighted residual methods used in the formulation of the finite element equations.