Solves the generalized Feynman-Kac PDE on a rectangular grid using a finite element Galerkin method. Initial and boundary conditions are provided. The solution is represented by a series of C2 Hermite quintic splines.

The generalized Feynman-Kac differential equation has the form

where the initial data satisfies

FEYNMAN_KAC uses a finite element Galerkin method over the rectangle

On each subinterval the solution is represented by

by

The Galerkin principle is then applied. Using the provided initial and boundary conditions leads to an index 1 differential-algebraic equation (DAE) for the time-dependent coefficients

This system is integrated using the variable order, variable step algorithm DASPH. Solution values and their time derivatives are returned at a grid preceding time T, expressed in units of time remaining.

More mathematical details are found in Hanson, R. (2008) “Integrating Feynman-Kac Equations Using Hermite Quintic Finite Elements”.

In Beckers (1980) there is a model for a Stochastic Differential Equation of option pricing. The idea is a “constant elasticity of variance diffusion (or CEV) class”

The Black-Scholes model is the limiting case α → 2. A numerical solution of this diffusion model yields the price of a call option. Various values of the strike price K, time values σ, and power coefficient α are used to evaluate the option price at values of the underlying price. The sets of parameters in the computation are:

With this model the Feynman-Kac differential equation is defined by identifying:

(Example feynman_kac_ex1.f90)

The value of the American Option on a Vanilla Put can be no smaller than its European counterpart. That is due to the American Option providing the opportunity to exercise at any time prior to expiration. This example compares this difference – or premium value of the American Option – at two time values using the Black-Scholes model. The example is based on Wilmott et al. (1996, p. 176), and uses the non-linear forcing or weighting term described in Hanson, R. (2008), “Integrating Feynman-Kac Equations Using Hermite Quintic Finite Elements”, for evaluating the price of the American Option. A call to the subroutine fkinit_put sets the initial conditions. One breakpoint is set exactly at the strike price.

The sets of parameters in the computation are:

(Example feynman_kac_ex2.f90)

This example evaluates the price of a European Option using two payoff strategies: Cash-or-Nothing and Vertical Spread. In the first case the payoff function is

The value B is regarded as the bet on the asset price, see Wilmott et al. (1995, p. 39-40). The second case has the payoff function

Both problems use the same boundary conditions. Each case requires a separate integration of the Black-Scholes differential equation, but only the payoff function evaluation differs in each case. The sets of parameters in the computation are:

(Example feynman_kac_ex3.f90)

This example evaluates the price of a convertible bond. Here, convertibility means that the bond may, at any time of the holder’s choosing, be converted to a multiple of the specified asset. Thus a convertible bond with price x returns an amount K at time T unless the owner has converted the bond to νx, ν ≥ 1, units of the asset at some time prior to T. This definition, the differential equation and boundary conditions are given in Chapter 18 of Wilmott et al. (1996). Using a constant interest rate and volatility factor, the parameters and boundary conditions are:

Note that the error tolerance is set to a pure absolute error of value 10-3. The free boundary constraint f(x,t) = νx is achieved by use of a non-linear forcing term in the subroutine fkforce_cbond. The terminal conditions are provided with the user subroutine fkinit_cbond.

(Example feynman_kac_ex4.f90)

This example illustrates a method for evaluating a certain “Bermudan Style” or non-standard American option. These options are American Style options restricted to certain dates where the option may be exercised. Since this agreement gives the holder more opportunity than a European option, it is worth more. But since the holder can only exercise at certain times it is worth no more than the American style option value that can be exercised at any time. Our solution method uses the same model and data as in Example 2, but allows exercise at weekly intervals. Thus we integrate, for half a year, over each weekly interval using a European style Black-Scholes model, but with initial data at each new week taken from the corresponding values of the American style option.

(Example feynman_kac_ex5.f90)

Our previous examples are from the field of financial engineering. A final example is a physical model. The Oxygen Diffusion Problem is summarized in Crank [4], p. 10-20, 261-262. We present the numerical treatment of the transformed one-dimensional system

A slight difference from the Crank development is that we have reflected the time variable t → -t to match our form of the Feynman-Kac equation. We have a free boundary problem because the interface s(t) is implicit. This interface is implicitly defined by solving the variational relation (ft + fxx ‑1)f = 0, f ≥ 0. The first factor is zero for 0 ≤ x < s(t) and the second factor is zero for s(t) ≤ x ≤ 1. We list the Feynman-Kac equation coefficients, forcing term and boundary conditions, followed by comments.

In the example code, values of f(0,t) are checked against published values for certain values of t. Also checked are values of f(0,s(t)) = 0 at published values of the free boundary, for the same values of t.

(Example feynman_kac_ex6.f90)

This example calculates FK and BS solutions V(S,t) to the BS problem and the following Greeks:

Intrinsic Greeks include derivatives involving only S and t, the intrinsic FK arguments. In the above list, Value, Delta, Gamma, Theta, Charm, Color and Speed are all intrinsic Greeks. As is discussed in Hanson, R. (2008) “Integrating Feynman-Kac Equations Using Hermite Quintic Finite Elements”, the expansion of the FK solution function V(S,t) in terms of quintic polynomial functions defined on S-grid subintervals and subject to continuity constraints in derivatives 0, 1 and 2 across the boundaries of these subintervals allows Value, Delta, Gamma, Theta, Charm and Color to be calculated directly by routines FEYNMAN_KAC and HQSVAL.

Non-intrinsic Greeks are derivatives of V involving FK parameters other than the intrinsic arguments S and t, such as r and α. Non-intrinsic Greeks in the above list include Rho, Vega, Volga and Vanna. In order to calculate non-intrinsic Greek (parameter) derivatives or intrinsic Greek S-derivatives beyond the second (such as Speed) or t-derivatives beyond the first, the entire FK solution must be calculated 3 times (for a parabolic fit) or five times (for a quartic fit), at the point where the derivative is to be evaluated and at nearby points located symmetrically on either side.

Using a Taylor series expansion of f(σ + ɛ) truncated to m + 1 terms (to allow an m-degree polynomial fit of m+1 data points),

we are able to derive the following parabolic (3 point) estimation of first and second derivatives f (1)(σ) and f (2)(σ) in terms of the three values f(σ ‑ ɛ ), f(σ) and f(σ + ɛ), where ɛ = ɛfracσ and 0 < ɛfrac << 1:

Similarly, the quartic (5 point) estimation of f (1)(σ) and f (2)(σ) in terms of f(σ ‑ 2ɛ), f(σ ‑ ɛ ), f(σ), f(σ + ɛ), and f(σ  + 2ɛ) is:

For our example, the quartic estimate does not appear to be significantly better than the parabolic estimate, so we have produced only parabolic estimates by setting variable iquart to 0. The user may try the example with the quartic estimate simply by setting iquart to 1.

As is pointed out in “Integrating Feynman-Kac Equations Using Hermite Quintic Finite Elements”, the quintic polynomial expansion function used by FEYNMAN_KAC only allows for continuous derivatives through the second derivative. While up to fifth derivatives can be calculated from the quintic expansion (indeed function HQSVAL will allow the third derivative to be calculated by setting optional argument IDERIV to 3, as is done in this example), the accuracy is compromised by the inherent lack of continuity across grid points (i.e. subinterval boundaries).

This procedure is implemented in the current example to calculate the Greek Speed. (For comparison purposes, Speed is also calculated directly by setting the optional argument IDERIV to 3 in the call to function HQSVAL. The output from this direct calculation is called “Speed2”.)

To reach better accuracy results, all computations are done in double precision.

The average and maximum relative errors (defined as the absolute value of the difference between the BS and FK values divided by the BS value) for each of the Greeks is given at the end of the output. (These relative error statistics are given for nine combinations of Strike Price and volatility, but only one of the nine combinations is actually printed in the output.) Both intrinsic and non-intrinsic Greeks have good accuracy (average relative error is in the range 0.01 – 0.00001) except for Volga, which has an average relative error of about 0.05. This is probably a result of the fact that Volga involves differences of differences, which will more quickly erode accuracy than calculations using only one difference to approximate a derivative. Possible ways to improve upon the 2 to 4 significant digits of accuracy achieved in this example include increasing FK integration accuracy by reducing the initial stepsize (via optional argument RINITSTEPSIZE), by choosing more closely spaced S and t grid points (by adjusting FEYNMAN_KAC’s input arrays XGRID and TGRID) and by adjusting ɛfrac so that the central differences used to calculate the derivatives are not too small to compromise accuracy.

(Example feynman_kac_ex7.f90)