feynmanKac

This function solves the generalized Feynman-Kac PDE on a rectangular grid using a finite element Galerkin method. Initial and boundary conditions are satisfied. The solution is represented by a series of \(C^2\) Hermite quintic splines.

Synopsis

feynmanKac (nlbcd, nrbcd, xgrid, tgrid, fcnFkcfiv, fcnFkbcp, y, yPrime)

Required Arguments

int nlbcd (Input)

Number of left boundary conditions. It is required that 1 ≤ nlbcd ≤ 3.

int nrbcd (Input)

Number of right boundary conditions. It is required that 1 ≤ nrbcd ≤ 3.

float xgrid[] (Input)

Array of length nxgrid containing the breakpoints for the Hermite quintic splines used in the x discretization. The points in xgrid must be strictly increasing. The values xgrid[0] and xgrid[nxgrid-1] are the endpoints of the interval.

float tgrid[] (Input)

Array of length ntgrid containing the set of time points (in time-remaining units) where an approximate solution is returned. The points in tgrid must be positive and strictly increasing.

void fcnFkcfiv (x, t, iflag, value)

User-supplied function to compute the values of the coefficients \(\sigma, \sigma', \mu, \kappa\) for the Feynman-Kac PDE and the initial data function \(p(x), x_{\min} \leq x \leq x_{\max}\) .

float x (Input)

Space variable.

float t (Input)

Time variable.

int iflag (Input/Output)

On entry, iflag indicates which coefficient or data function is to be computed. The following table shows which value has to be returned by fcnFkcfiv for all possible values of iflag:

iflag	Computed coefficient
-1	\(\sigma' = \tfrac{\partial \sigma(x,t)}{\partial x}\)
0	\(p(x)\)
1	\(\sigma\)
2	\(\mu\)
3	\(\kappa\)

For non-zero input values of iflag, note when a coefficient does not depend on t. This is done by setting iflag = 0 after the coefficient is defined. If there is time dependence, the value of iflag should not be changed. This action will usually yield a more efficient algorithm because some finite element matrices do not have to be reassembled for each t value.

float value (Output)

Value of the coefficient or initial data function. Which value is computed depends on the input value for iflag, see description of iflag.

void fcnFkbcp (nbc, t, iflag, values[])

User-supplied function to define boundary values that the solution of the differential equation must satisfy. There are nlbcd conditions specified at the left end, \(x_{\text{min}}\), and nrbcd conditions at the right end, \(x_{\text{max}}\). The boundary conditions are

\[a(x,t)f + b(x,t)f_x + c(x,t)f_{xx} = d(x,t), \phantom{...} x = x_{\min} \text{ or } x = x_{\max}\]

int nbc (Input)

Number of boundary conditions. At \(x_{min}\), nbc=nlbcd, at \(x_{max}\), nbc = nrbcd.

float t (Input)

Time point of the boundary conditions.

int iflag (Input/Output)

On entry, iflag indicates whether the coefficients for the left or right boundary conditions are to be computed:

iflag	Computed boundary conditions
1	Left end, \(x = x_{min}\)
0	Right end, \(x = x_{max}\)

If there is no time dependence for one of the boundaries then set iflag = 0 after the array values is defined for either end point. This flag can avoid unneeded continued computation of the finite element matrices.

float values[] (Output)

Array of length 4 * max (nlbcd, nrbcd) containing the values of the boundary condition coefficients in its first 4*nbc locations. The coefficients for \(x_{min}\)are stored row-wise according to the following scheme:

\[\begin{split}\left( \begin{array}{c} a_1\left(x_{\min}, t\right), b_1\left(x_{\min}, t\right), c_1\left(x_{\min}, t\right), d_1\left(x_{\min}, t\right) \\ \vdots \\ a_{\mathit{nlbcd}}\left(x_{\min}, t\right), b_{\mathit{nlbcd}}\left(x_{\min}, t\right), c_{\mathit{nlbcd}}\left(x_{\min}, t\right), d_{\mathit{nlbcd}}\left(x_{\min}, t\right) \\ \end{array} \right)\end{split}\]

The coefficients for \(x_{max}\) are stored similarly.

float y[] (Output)

An array of size (ntgrid+1) by (3 × nxgrid) containing the coefficients of the Hermite representation of the approximate solution for the Feynman-Kac PDE at time points 0, tgrid[0], …, tgrid[ntgrid-1]. The approximate solution is given by

\[f(x,t) = \sum_{j=0}^{3 * \mathtt{nxgrid}-1} y_{ij} \beta_j(x)\]

for

\[t = \mathtt{tgrid}[i-1], i = 1, \ldots, \mathtt{ntgrid}\]

The representation for the initial data at \(t = 0\) is

\[p(x) = \sum_{j=0}^{3 * \mathtt{nxgrid}-1} y_{0j} \beta_j(x)\]

The (ntgrid + 1) by (3 * nxgrid) matrix

\[\left(y_{ij}\right)_{i=0,\ldots,\mathtt{ntgrid}}^{j=0,\ldots,3*\mathtt{nxgrid}-1}\]

is stored row-wise in array y.

After the integration, use row i of array y as input argument coef in function feynmanKacEvaluate to obtain an array of values for f(x, t) or its partials \(f_x, f_{xx}, f_{xxx}\) at time point t=0, tgrid[i-1], i=1,…, ntgrid.

The expressions for the basis functions \(\beta_j(x)\) are represented piece-wise and can be found in Hanson, R. (2008) “Integrating Feynman-Kac Equations Using Hermite Quintic Finite Elements.”

float yPrime[] (Output)

An array of size (ntgrid + 1) by (3 × nxgrid) containing the first derivatives of y at time points 0, tgrid[0], ..., tgrid[ntgrid - 1], i.e.,

\[f_t\left(x,\overline{t}\right) = \sum_{j=0}^{3*\mathtt{nxgrid}-1} y'_{ij} \beta_j(x)\]

for

\[\overline{t} = \mathtt{tgrid}[i-1], i=1, \ldots, \mathtt{ntgrid}\]

and

\[f_t\left(x,\overline{t}\right) = \sum_{j=0}^{3*\mathtt{nxgrid}-1} y'_{0j} \beta_j(x) \phantom{...} \text{ for } \overline{t} = 0 \text{.}\]

The (ntgrid + 1) by (3 \*nxgrid) matrix

\[\left(y'_{ij}\right)_{i=0,\ldots,\mathtt{ntgrid}}^{j=0,\ldots,3*\mathtt{nxgrid}-1}\]

is stored row-wise in array yPrime.

After the integration, use row i of array yPrime as input argument coef in function feynmanKacEvaluate to obtain an array of values for the partials \(f_t, f_{tx}, f_{txx}, f_{txxx}\) at time point t=tgrid[i-1], i=1, ..., ntgrid, and row 0 for the partials at t=0.

Optional Arguments

fcnInit, void fcnInit(nxgrid, ntgrid, xgrid[], tgrid[], time, yprime[], y[], atol[], rtol[]) (Input)

User-supplied function for adjustment of initial data or as an opportunity for output during the integration steps. The solution values of the model parameters are presented in the arrays y and yprime at the integration points time = 0, tgrid[j], j=0,…, ntgrid-1. At the initial point, integral least-squares estimates are made for representing the initial data \(p(x)\). If this is not satisfactory, fcnInit can change the contents of y[] to match data for any reason.

int nxgrid (Input)

Number of grid lines in the x-direction.

int ntgrid (Input)

Number of time points where an approximate solution is returned.

float xgrid[] (Input)

Vector of length nxgrid containing the breakpoints for the Hermite quintic splines used in the x discretization.

float tgrid[] (Input)

Vector of length ntgrid containing the set of time points (in time-remaining units) where an approximate solution is returned.

float time (Input)

Time variable.

float yprime[] (Input)

Vector of length 3\*nxgrid containing the derivative of the coefficients of the Hermite quintic spline at time point time.

float y[] (Input/Output)

Vector of length 3 × nxgrid containing the coefficients of the Hermite quintic spline at time point time.

float atol[] (Input/Output)

Array of length 3 × nxgridcontaining absolute error tolerances.

float rtol[] (Input/Output)

Array of length 3 × nxgrid containing relative error tolerances.

fcnForce, void fcnForce(interval, ndeg, nxgrid, y[], time, width, xlocal[], qw[], u[], phi[], dphi[]) (Input)

Function fcnForce is used in case there is a non-zero term \(\phi(f, x, t)\) in the Feynman-Kac differential equation. If function fcnForce is not used, it is assumed that \(\phi(f, x, t)\) is identically zero.

int interval (Input)

Index indicating the bounds xgrid[interval-1] and xgrid[interval] of the integration interval, 1≤interval ≤ nxgrid-1.

int ndeg (Input)

The degree used for the Gauss-Legendre formulas.

int nxgrid (Input)

Number of grid lines in the x-direction.

float y[] (Input)

Vector of length 3\*nxgrid containing the coefficients of the Hermite quintic spline representing the solution of the Feynman-Kac equation at time point time.

float time (Input)

Time variable.

float width (Input)

The interval width, width = xgrid[interval] - xgrid[interval-1].

float xlocal[] (Input)

Array of length ndeg containing the Gauss-Legendre points translated and normalized to the interval [xgrid[interval-1], xgrid[interval]].

float qw[] (Input)

Vector of length ndeg containing the Gauss-Legendre weights.

float u[] (Input)

Array of dimension 12 by ndeg containing the basis function values that define \(\tilde{\beta}(x)\) at the Gauss-Legendre points xlocal. Setting \(u_{k,i} :=\) u[i+k\*ndeg] and \(x_i :=\) xlocal[i], \(\tilde{\beta}\left( x_i \right)\) is defined as

\[\begin{split}\begin{array}{l} \tilde{\beta}\left(x_i\right): = \left(\beta_{3 * (\text{interval} - 1)}\left(x_i\right), \ldots, \beta_{3 * (\text{interval} + 2)}\left(x_i\right)\right)^T \\ = \left(u_{0,i},u_{1,i},u_{2,i},u_{6,i},u_{7,i},u_{8,i}\right)^T. \end{array}\end{split}\]

float phi[] (Output)

Vector of length 6 containing Gauss-Legendre approximations for the local contributions

\[\varphi_{\mathrm{t}} := \int_{\mathrm{xgrid[interval-1]}}^{\mathrm{xgrid[interval]}} \mathit{\Phi}(f,x,t) \tilde{\beta}(x) dx,\]

where \(t =\) time and

\[\tilde{\beta}(x) := \left(\beta_{3*(\mathrm{interval}-1)}(x), \ldots, \beta_{3*(\mathrm{interval}+2)}(x)\right)^T\]

Vector phi contains elements

\[\mathrm{phi}[i] = \mathrm{width} * \textstyle\sum_{j=0}^{\mathrm{ndegGL}-1} qw[j] \tilde{\beta}_i \left(x_j\right) \phi(f, \mathrm{xlocal}[j], \mathrm{time})\]

for i=0,…,5.

float dphi[] (Output)

Array of dimension 6 by 6 containing a Gauss-Legendre approximation for the Jacobian of the local contributions \(\phi_t\) at \(t =\) time,

\[\frac{\partial \varphi_t}{\partial y} = \int_{\mathrm{xgrid[interval-1]}}^{\mathrm{xgrid[interval]}} \frac{\partial \mathit{\Phi}(f,x,t)}{\partial f} \tilde{\beta}(x)\tilde{\beta}^T(x) dx.\]

The approximation to this symmetric matrix is stored row-wise, i.e.

\[\mathrm{dphi}[j+i*s] = \mathrm{width} * \textstyle\sum_{k=0}^{n\mathrm{deg}GL-1} qw[k]\phantom{...} \tilde{\beta}_i\left(x_k\right) \tilde{\beta}_j\left(x_{\mathrm{k}}\right) \frac{\partial \mathit{\Phi}}{\partial f} |_{x=\mathit{xlocal}[\mathrm{k}],t=\mathit{time}}\]

for \(i, j = 0, …, 5\).

atolRtolScalars, float atol, float rtol (Input)

Scalar values that apply to the error estimates of all components of the solution y in the differential equation solver SDASLX. See optional argument atolRtolArrays if separate tolerances are to be applied to each component of y.

Default: atol and rtol are set to \(10^{-3}\) in single precision and \(10^{-5}\) in double precision.

atolRtolArrays, float atol[], float rtol[], (Input)

Componentwise tolerances are used for the computation of solution y in the differential equation solver SDASLX. Arguments atol and rtol are arrays of length 3 × nxgrid to be used for the absolute and relative tolerance and to be applied to each component of the solution, y. See optional argument atolRtolScalars if scalar values of absolute and relative tolerances are to be applied to all components.

Default: All elements of atol and rtol are set to \(10^{-3}\) in single precision and \(10^{-5}\) in double precision.

nDegree, int (Input)

The degree used for the Gauss-Legendre formulas for constructing the finite element matrices. It is required that nDegree ≥ 6.

Default: nDegree = 6.

tdepend (Output)

Vector of length 7 indicating time dependence of the coefficients, boundary conditions and function φ in the Feynman-Kac equation. If tdepend[i] = 0 then argument i is not time dependent, if tdepend[i]=1 then argument i is time dependent.

i	Computed coefficient
0	\(\sigma'\)
1	\(\sigma\)
2	\(\mu\)
3	\(\kappa\)
4	Left boundary conditions
5	Right boundary conditions
6	φ

maxStep, float (Input)

This is the maximum step size the integrator may take.

Default: maxStep = machine(2), the largest possible machine number.

initialStepsize, float (Input)

The starting step size for the integration. Must be less than zero since the integration is internally done from t=0 to t=tgrid[ntgrid-1] in a negative direction.

Default: initialStepsize = -Ɛ, where Ɛ is the machine precision.

maxNumberSteps, int (Input)

The maximum number of steps between each output point of the integration.

Default: maxNumberSteps = 500000.

stepControl, int (Input)

Indicates which step control algorithm is used for the integration. If stepControl = 0, then the step control method of Söderlind is used. If stepControl = 1, then the method used by the original Petzold code SASSL is used.

Default: stepControl = 0.

maxBdfOrder, int (Input)

Maximum order of the backward differentiation formulas used in the integrator. It is required that 1≤ maxBdfOrder ≤ 5.

Default: maxBdfOrder = 5.

tBarrier, float (Input)

This optional argument controls whether the code should integrate past a special point, tBarrier, and then interpolate to get \(y\) and \(y'\) at the points in tgrid[]. If this optional argument is present, the integrator assumes the equations either change on the alternate sides of tBarrier or they are undefined there. In this case, the code creeps up to tBarrier in the direction of integration. It is required that tBarrier ≥tgrid[ntgrid-1].

Default: tBarrier = tgrid[ntgrid-1].

istate (Output)

An array of length 7 whose entries flag the state of computation for the matrices and vectors required in the integration. For each entry, a zero indicates that no computation has been done or there is a time dependence. A one indicates that the entry has been computed and there is no time dependence. The istate entries are as follows:

i	Entry in istate
0	Mass Matrix, M
1	Stiffness matrix, N
2	Bending matrix, R
3	Weighted mass, K
4	Left boundary conditions
5	Right- boundary conditions
6	Forcing term

Default: istate[i] = 0 for I = 0,…,6.

evals (Output)

Array of length 3 summarizing the number of evaluations required during the integration.

i	evals[i]
0	Number of residual function evaluations of the DAE used in the model.
1	Number of factorizations of the differential matrix associated with solving the DAE.
2	Number of linear system solve steps using the differential matrix.

Description

The generalized Feynman-Kac differential equation has the form

\[f_t + \mu(x,t)f_x + \frac{\sigma^2(x,t)}{2}f_{xx} - \kappa(x,t)f = \phi(f,x,t),\]

where the initial data satisfies \(f(x,T) = p(x)\). The derivatives are

\[f_t = \frac{\partial f}{\partial t}\]

The function feynmanKac uses a finite element Galerkin method over the rectangle

\[\left[x_{\min},x_{\max}\right] \times \left[\hat{T},T\right]\]

in \((x,t)\) to compute the approximate solution. The interval \(\left[ x_{\min}, x_{\max} \right]\) is decomposed with a grid

\[\left(x_{\min} =\right) x_1 < x_2 < \ldots < x_m \left(= x_{\max}\right).\]

On each subinterval the solution is represented by

\[\begin{split}\begin{array}{c} f(x,t) = f_ib_0(z) + f_{i+1}b_0(1-z) + h_if'_ib_1(z) - h_if'_{i+1}b_1(1-z) \\ + h_i^2f''_ib_2(z) + h_i^2f''_i+1b_2(1-z). \end{array}\end{split}\]

The values

\[f_i,f'_i,f''_i,f_{i+1},f'_{i+1},f''_{i+1}\]

are time-dependent coefficients associated with each interval. The basis functions \(b_0, b_1, b_2\) are given for

\[x \in \left[x_i,x_{i+1}\right], h_i = x_{i+1} - x_i, z = \left(x-x_i\right) / h_i, z \in [0,1],\]

by

\[\begin{split}\begin{aligned} b_0(z) &= -6z^5 + 15z^4 - 10z^3 + 1 = (1-z)^3\left(6z^2 + 3z + 1\right), \\ b_1(z) &= -3z^5 + 8z^4 - 6z^3 + z = (1-z)^3z(3z+1), \\ b_2(z) &= \tfrac{1}{2}\left(-z^5+3z^4-3z^3+z^2\right) = \tfrac{1}{2}(1-z)^3z^2. \end{aligned}\end{split}\]

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

\[f_i,f'_i,f''_i,f_{i+1},f'_{i+1},f''_{i+1} \left[ = y_i,y_{i+1},y_{i+2}, \dots \right], i=1, \dots, m-1\]

This system is integrated using the variable order, variable step algorithm DDASLX/SDASLX noted in Hanson and Krogh, R. (2008) Solving Constrained Differential-Algebraic Systems Using Projections. Solution values and their time derivatives are returned at a grid preceding time T, expressed in units of time remaining. For further details of deriving and solving the system see Hanson, R. (2008) Integrating Feynman-Kac Equations Using Hermite Quintic Finite Elements. To evaluate f or its partial derivatives at any time point in the grid, use the function feynmanKacEvaluate.

Examples

Example 1

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. The coefficients, payoff, boundary conditions and forcing term for American and European options are defined by the user functions fkcfiv_put, fkbcp_put and fkforce_put, respectively. One breakpoint is set exactly at the strike price.

The sets of parameters in the computation are:

1. Strike price \(K = \{ 10.0 \}\).
2. Volatility \(\sigma = \{0.4\}\).
3. Times until expiration \(= \{ 1/4, 1/2 \}\).
4. Interest rate \(r = 0.1\).
5. \(x_{\min} = 0.0, \phantom{...} x_{\max} = 30.0\).
6. \(nx = 61, \phantom{...} n = 3 \times nx = 183\).

The payoff function is the “vanilla option”, \(p(x) = \max (K - x, 0)\).

from __future__ import print_function
from numpy import empty, zeros
from pyimsl.math.feynmanKac import feynmanKac
from pyimsl.math.feynmanKacEvaluate import feynmanKacEvaluate

# These routines define the coefficients, payoff, boundary conditions
# and forcing term for American and European Options.


def fkcfiv_put(x, tx, iflag, value):
    if iflag[0] == 0:
        # The payoff function
        value[0] = max(strike_price - x, 0.0)
    elif iflag[0] == -1:
        # The coefficient derivative d sigma/ dx
        value[0] = sigma
    elif iflag[0] == 1:
        # The coefficient sigma(x)
        value[0] = sigma * x
    elif iflag[0] == 2:
        # The coefficient mu(x)
        value[0] = (interest_rate - dividend) * x
    elif iflag[0] == 3:
        # The coefficient kappa(x)
        value[0] = interest_rate

    # Note that there is no time dependence
    iflag[0] = 0


def fkbcp_put(nbc, tx, iflag, val):
    if iflag[0] == 1:
        val[0] = 0.0
        val[1] = 1.0
        val[2] = 0.0
        val[3] = -1.0
        val[4] = 0.0
        val[5] = 0.0
        val[6] = 1.0
        val[7] = 0.0
    elif iflag[0] == 2:
        val[0] = 1.0
        val[1] = 0.0
        val[2] = 0.0
        val[3] = 0.0
        val[4] = 0.0
        val[5] = 1.0
        val[6] = 0.0
        val[7] = 0.0
        val[8] = 0.0
        val[9] = 0.0
        val[10] = 1.0
        val[11] = 0.0
    # Note no time dependence
    iflag[0] = 0


def fkforce_put(interval, ndeg, nxgrid, y, time,
                width, xlocal, qw, u, phi, dphi):
    local = 6
    yl = empty(6)
    bf = empty(6)
    for i in range(0, local):
        yl[i] = y[3 * interval - 3 + i]
        phi[i] = 0.0
    mu = 2
    # This is the local definition of the forcing term
    for j in range(1, local + 1):
        for l in range(1, ndeg + 1):
            bf[0] = u[(l - 1)]
            bf[1] = u[(l - 1) + ndeg]
            bf[2] = u[(l - 1) + 2 * ndeg]
            bf[3] = u[(l - 1) + 6 * ndeg]
            bf[4] = u[(l - 1) + 7 * ndeg]
            bf[5] = u[(l - 1) + 8 * ndeg]

            rt = 0.0
            for k in range(0, local):
                rt += yl[k] * bf[k]

            rt = value / (rt + value - (strike_price - xlocal[l - 1]))
            phi[j - 1] += qw[l - 1] * bf[j - 1] * pow(rt, mu)

    for i in range(0, local):
        phi[i] = -phi[i] * width * interest_rate * strike_price

    # This is the local derivative matrix for the forcing term
    for i in range(0, local * local):
        dphi[i] = 0.0

    for j in range(1, local + 1):
        for i in range(1, local + 1):
            for l in range(1, ndeg + 1):
                bf[0] = u[(l - 1)]
                bf[1] = u[(l - 1) + ndeg]
                bf[2] = u[(l - 1) + 2 * ndeg]
                bf[3] = u[(l - 1) + 6 * ndeg]
                bf[4] = u[(l - 1) + 7 * ndeg]
                bf[5] = u[(l - 1) + 8 * ndeg]

                rt = 0.0
                for k in range(0, local):
                    rt += yl[k] * bf[k]

                rt = 1.0 / (rt + value - (strike_price - xlocal[l - 1]))
                dphi[i - 1 + (j - 1) * local] += qw[l - 1] * \
                    bf[i - 1] * bf[j - 1] * pow(rt, mu + 1)

    for i in range(0, local * local):
        dphi[i] = mu * dphi[i] * width * \
            pow(value, mu) * interest_rate * strike_price

# Compute American Option Premium for Vanilla Put


NXGRID = 61
NTGRID = 2
NV = 9

# The strike price
strike_price = 10.0

# The sigma value
sigma = 0.4

# Time values for the options
nt = 2
time = [0.25, 0.5]

# Values of the underlying where evaluations are made
xs = [0.0, 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0]

# Value of the interest rate and continuous dividend
interest_rate = 0.1
dividend = 0.0

# Values of the min and max underlying values modeled
x_min = 0.0
x_max = 30.0

# Define parameters for the integration step.
nint = NXGRID - 1
n = 3 * NXGRID

xgrid = empty(NXGRID)
ye = empty(0)
yeprime = empty(0)
ya = empty(0)
yaprime = empty(0)
fe = zeros([nt, NV])
fa = zeros([nt, NV])

nlbcd = 2
nrbcd = 3
value = 0.2e-2

# Define an equally-spaced grid of points for the
# underlying price
dx = (x_max - x_min) / float(nint)
xgrid[0] = x_min
xgrid[NXGRID - 1] = x_max
for i in range(2, NXGRID):
    xgrid[i - 1] = xgrid[i - 2] + dx

feynmanKac(nlbcd, nrbcd, xgrid, time,
           fkcfiv_put, fkbcp_put, ye, yeprime,
           atolRtolScalars={'atol': 0.5e-2, 'rtol': 0.5e-2})

feynmanKac(nlbcd, nrbcd, xgrid, time,
           fkcfiv_put, fkbcp_put, ya, yaprime,
           fcnForce=fkforce_put,
           atolRtolScalars={'atol': 0.5e-2, 'rtol': 0.5e-2})

# Evaluate solutions at vector of points XS(:), at each
# time value prior to expiration.
for i in range(0, nt):
    tmp = feynmanKacEvaluate(xgrid, xs, ye[i + 1, ])
    fe[i, :] = tmp
    tmp = feynmanKacEvaluate(xgrid, xs, ya[i + 1, ])
    fa[i, :] = tmp

print("\nAmerican Option Premium for Vanilla Put, 3 and 6 Months Prior to Expiry")
print("       Number of equally spaced spline knots:%4d" % NXGRID)
print("       Number of unknowns:%4d" % n)
print("       Strike=%6.2f, sigma=%5.2f, Interest Rate=%5.2f\n" %
      (strike_price, sigma, interest_rate))
print("       Underlying            European            American")
for i in range(0, NV):
    print("       %10.4f%10.4f%10.4f%10.4f%10.4f" % (xs[i], fe[0, i],
                                                     fe[1, i], fa[0, i], fa[1, i]))

Output

American Option Premium for Vanilla Put, 3 and 6 Months Prior to Expiry
       Number of equally spaced spline knots:  61
       Number of unknowns: 183
       Strike= 10.00, sigma= 0.40, Interest Rate= 0.10

       Underlying            European            American
           0.0000    9.7535    9.5139   10.0000   10.0000
           2.0000    7.7535    7.5139    8.0000    8.0000
           4.0000    5.7536    5.5159    6.0000    6.0000
           6.0000    3.7607    3.5691    4.0000    4.0000
           8.0000    1.9061    1.9168    2.0174    2.0856
          10.0000    0.6560    0.8533    0.6780    0.9035
          12.0000    0.1638    0.3362    0.1681    0.3521
          14.0000    0.0364    0.1269    0.0371    0.1322
          16.0000    0.0083    0.0486    0.0083    0.0501

Example 2

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”

\[dS = \mu S dt + \sigma S^{\alpha/2} dW, \phantom{...} 0 \leq \alpha < 2\]

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

1. power \(\alpha = \{2.0, 1.0, 0.0\}\).
2. strike price \(K = \{15.0, 20.0, 25.0\}\).
3. volatility \(\sigma = \{0.2, 0.3, 0.4\}\).
4. times until expiration \(= \{ 1/12, 4/12, 7/12 \}\).
5. underlying prices \(= \{ 19.0, 20.0, 21.0 \}\).
6. interest rate \(r = 0.05\).
7. \(x_{\min} = 0, x_{\max} = 60\).
8. \(nx = 121, \phantom{...} n = 3 \times nx = 363\).

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

\[x: S\]

\[\sigma(x,t): \phantom{...} \sigma x^{\alpha/2}; \phantom{...} \frac{\partial \sigma}{\partial x} = \frac{a\sigma}{2} x^{\alpha/2-1}\]

\[\mu(x,t): \phantom{...} rx\]

\[\kappa(x,t): \phantom{...} r\]

\[\phi(f,x,t) \equiv 0\]

The payoff function is the “vanilla option”, \(p(x) = \max (x - K, 0)\).

from __future__ import print_function
from numpy import *
from pyimsl.math.feynmanKac import feynmanKac
from pyimsl.math.feynmanKacEvaluate import feynmanKacEvaluate


def fcn_fkcfiv(x, tx, iflag, value):
    if iflag[0] == 0:
        # The payoff function
        value[0] = max(x - strike_price, 0.0)
    elif iflag[0] == -1:
        # The coefficient derivative d sigma/ dx
        value[0] = 0.5 * alpha * sigma * pow(x, alpha * 0.5 - 1.0)
    elif iflag[0] == 1:
        # The coefficient sigma(x)
        value[0] = sigma * pow(x, alpha * 0.5)
    elif iflag[0] == 2:
        # The coefficient mu(x)
        value[0] = (interest_rate - dividend) * x
    elif iflag[0] == 3:
        # The coefficient kappa(x)
        value[0] = interest_rate

    # Note that there is no time dependence
    iflag[0] = 0


def fcn_fkbcp(nbc, tx, iflag, val):
    if iflag[0] == 1:
        val[0] = 1.0
        val[1] = 0.0
        val[2] = 0.0
        val[3] = 0.0
        val[4] = 0.0
        val[5] = 1.0
        val[6] = 0.0
        val[7] = 0.0
        val[8] = 0.0
        val[9] = 0.0
        val[10] = 1.0
        val[11] = 0.0
        # Note no time dependence
        iflag[0] = 0
    elif iflag[0] == 2:
        df = exp(interest_rate * tx)
        val[0] = 1.0
        val[1] = 0.0
        val[2] = 0.0
        val[3] = x_max - df * strike_price
        val[4] = 0.0
        val[5] = 1.0
        val[6] = 0.0
        val[7] = 1.0
        val[8] = 0.0
        val[9] = 0.0
        val[10] = 1.0
        val[11] = 0.0


# Compute Constant Elasticity of Variance Model for Vanilla Call
NXGRID = 121
NTGRID = 3
NV = 3

# The set of strike prices
KS = [15.0, 20.0, 25.0]

# The set of sigma values
sigmas = [0.2, 0.3, 0.4]

# The set of model diffusion powers
alphas = [2.0, 1.0, 0.0]

# Time values for the options
nt = 3
time = [1.0 / 12.0, 4.0 / 12.0, 7.0 / 12.0]

# Values of the underlying where evaluations are made
xs = [19.0, 20.0, 21.0]

# Value of the interest rate and continuous dividend
interest_rate = 0.05
dividend = 0.0

# Values of the min and max underlying values modeled
x_min = 0.0
x_max = 60.0

# Define parameters for the integration step.
nint = NXGRID - 1
n = 3 * NXGRID
xgrid = empty(NXGRID)
y = empty(0)
yprime = empty(0)
f = empty([NTGRID, NV])

# Number of left/right boundary conditions
nlbcd = 3
nrbcd = 3

# Define equally-spaced grid of points for the underlying price
dx = (x_max - x_min) / float(nint)
xgrid[0] = x_min
xgrid[NXGRID - 1] = x_max

for i in range(2, NXGRID):
    xgrid[i - 1] = xgrid[i - 2] + dx

print("  Constant Elasticity of Variance Model for Vanilla Call")
print("       Interest Rate:%7.3f\tContinuous Dividend:%7.3f" %
      (interest_rate, dividend))
print("       Minimum and Maximum Prices of Underlying:%7.2f%7.2f" % (x_min, x_max))
print("       Number of equally spaced spline knots:%4d" % (NXGRID - 1))
print("       Number of unknowns:%4d\n" % n)
print("       Time in Years Prior to Expiration:%7.4f%7.4f%7.4f"
      % (time[0], time[1], time[2]))
print("       Option valued at Underlying Prices:%7.2f%7.2f%7.2f\n"
      % (xs[0], xs[1], xs[2]))

for i1 in range(1, 4):         # Loop over power
    for i2 in range(1, 4):     # Loop over volatility
        for i3 in range(1, 4):  # Loop over strike price
            # Pass data through into evaluation routines.
            strike_price = KS[i3 - 1]
            sigma = sigmas[i2 - 1]
            alpha = alphas[i1 - 1]

            feynmanKac(nlbcd, nrbcd, xgrid, time,
                       fcn_fkcfiv, fcn_fkbcp, y, yprime)

            # Evaluate solution at vector of points xs, at each time
            # Value prior to expiration.
            for i in range(0, NTGRID):
                tmp = feynmanKacEvaluate(xgrid, xs, y[i + 1, ])
                f[i, :] = tmp

            print("  Strike=%5.2f, Sigma=%5.2f, Alpha=%5.2f" %
                  (KS[i3 - 1], sigmas[i2 - 1], alphas[i1 - 1]))

            for i in range(0, NV):
                print("                       Call Option Values", end=' ')
                for j in range(0, nt):
                    print("%7.4f  " % f[j, i], end=' ')
                print()
            print()

Output

  Constant Elasticity of Variance Model for Vanilla Call
       Interest Rate:  0.050	Continuous Dividend:  0.000
       Minimum and Maximum Prices of Underlying:   0.00  60.00
       Number of equally spaced spline knots: 120
       Number of unknowns: 363

       Time in Years Prior to Expiration: 0.0833 0.3333 0.5833
       Option valued at Underlying Prices:  19.00  20.00  21.00

  Strike=15.00, Sigma= 0.20, Alpha= 2.00
                       Call Option Values  4.0624    4.2576    4.4734   
                       Call Option Values  5.0624    5.2505    5.4492   
                       Call Option Values  6.0624    6.2486    6.4386   

  Strike=20.00, Sigma= 0.20, Alpha= 2.00
                       Call Option Values  0.1312    0.5951    0.9693   
                       Call Option Values  0.5024    1.0880    1.5093   
                       Call Option Values  1.1980    1.7478    2.1745   

  Strike=25.00, Sigma= 0.20, Alpha= 2.00
                       Call Option Values  0.0000    0.0111    0.0751   
                       Call Option Values  0.0000    0.0376    0.1630   
                       Call Option Values  0.0006    0.1036    0.3150   

  Strike=15.00, Sigma= 0.30, Alpha= 2.00
                       Call Option Values  4.0636    4.3405    4.6627   
                       Call Option Values  5.0625    5.2951    5.5794   
                       Call Option Values  6.0624    6.2712    6.5248   

  Strike=20.00, Sigma= 0.30, Alpha= 2.00
                       Call Option Values  0.3107    1.0261    1.5479   
                       Call Option Values  0.7317    1.5404    2.0999   
                       Call Option Values  1.3762    2.1674    2.7362   

  Strike=25.00, Sigma= 0.30, Alpha= 2.00
                       Call Option Values  0.0005    0.1124    0.3564   
                       Call Option Values  0.0035    0.2184    0.5565   
                       Call Option Values  0.0184    0.3869    0.8230   

  Strike=15.00, Sigma= 0.40, Alpha= 2.00
                       Call Option Values  4.0755    4.5143    4.9673   
                       Call Option Values  5.0660    5.4210    5.8328   
                       Call Option Values  6.0633    6.3588    6.7306   

  Strike=20.00, Sigma= 0.40, Alpha= 2.00
                       Call Option Values  0.5109    1.4625    2.1260   
                       Call Option Values  0.9611    1.9934    2.6915   
                       Call Option Values  1.5807    2.6088    3.3202   

  Strike=25.00, Sigma= 0.40, Alpha= 2.00
                       Call Option Values  0.0081    0.3302    0.7795   
                       Call Option Values  0.0287    0.5178    1.0656   
                       Call Option Values  0.0820    0.7690    1.4097   

  Strike=15.00, Sigma= 0.20, Alpha= 1.00
                       Call Option Values  4.0624    4.2479    4.4312   
                       Call Option Values  5.0624    5.2479    5.4312   
                       Call Option Values  6.0624    6.2479    6.4312   

  Strike=20.00, Sigma= 0.20, Alpha= 1.00
                       Call Option Values  0.0000    0.0219    0.1051   
                       Call Option Values  0.1497    0.4107    0.6484   
                       Call Option Values  1.0832    1.3314    1.5773   

  Strike=25.00, Sigma= 0.20, Alpha= 1.00
                       Call Option Values -0.0000   -0.0000    0.0000   
                       Call Option Values -0.0000   -0.0000    0.0000   
                       Call Option Values -0.0000   -0.0000    0.0000   

  Strike=15.00, Sigma= 0.30, Alpha= 1.00
                       Call Option Values  4.0624    4.2479    4.4312   
                       Call Option Values  5.0624    5.2479    5.4312   
                       Call Option Values  6.0624    6.2479    6.4312   

  Strike=20.00, Sigma= 0.30, Alpha= 1.00
                       Call Option Values  0.0010    0.0786    0.2208   
                       Call Option Values  0.1993    0.4997    0.7539   
                       Call Option Values  1.0835    1.3444    1.6022   

  Strike=25.00, Sigma= 0.30, Alpha= 1.00
                       Call Option Values -0.0000    0.0000    0.0000   
                       Call Option Values -0.0000    0.0000    0.0000   
                       Call Option Values -0.0000    0.0000    0.0004   

  Strike=15.00, Sigma= 0.40, Alpha= 1.00
                       Call Option Values  4.0624    4.2479    4.4312   
                       Call Option Values  5.0624    5.2479    5.4312   
                       Call Option Values  6.0624    6.2479    6.4312   

  Strike=20.00, Sigma= 0.40, Alpha= 1.00
                       Call Option Values  0.0072    0.1540    0.3446   
                       Call Option Values  0.2498    0.5950    0.8728   
                       Call Option Values  1.0868    1.3795    1.6586   

  Strike=25.00, Sigma= 0.40, Alpha= 1.00
                       Call Option Values  0.0000    0.0000    0.0000   
                       Call Option Values  0.0000    0.0000    0.0005   
                       Call Option Values  0.0000    0.0002    0.0057   

  Strike=15.00, Sigma= 0.20, Alpha= 0.00
                       Call Option Values  4.0624    4.2479    4.4311   
                       Call Option Values  5.0624    5.2479    5.4312   
                       Call Option Values  6.0624    6.2479    6.4312   

  Strike=20.00, Sigma= 0.20, Alpha= 0.00
                       Call Option Values  0.0001    0.0001    0.0002   
                       Call Option Values  0.0816    0.3316    0.5748   
                       Call Option Values  1.0817    1.3308    1.5748   

  Strike=25.00, Sigma= 0.20, Alpha= 0.00
                       Call Option Values  0.0000   -0.0000   -0.0000   
                       Call Option Values  0.0000   -0.0000   -0.0000   
                       Call Option Values -0.0000    0.0000   -0.0000   

  Strike=15.00, Sigma= 0.30, Alpha= 0.00
                       Call Option Values  4.0624    4.2479    4.4312   
                       Call Option Values  5.0624    5.2479    5.4312   
                       Call Option Values  6.0624    6.2479    6.4312   

  Strike=20.00, Sigma= 0.30, Alpha= 0.00
                       Call Option Values  0.0000   -0.0000    0.0026   
                       Call Option Values  0.0894    0.3326    0.5753   
                       Call Option Values  1.0826    1.3306    1.5749   

  Strike=25.00, Sigma= 0.30, Alpha= 0.00
                       Call Option Values  0.0000   -0.0000    0.0000   
                       Call Option Values  0.0000   -0.0000    0.0000   
                       Call Option Values  0.0000   -0.0000   -0.0000   

  Strike=15.00, Sigma= 0.40, Alpha= 0.00
                       Call Option Values  4.0624    4.2479    4.4312   
                       Call Option Values  5.0624    5.2479    5.4312   
                       Call Option Values  6.0624    6.2479    6.4312   

  Strike=20.00, Sigma= 0.40, Alpha= 0.00
                       Call Option Values -0.0000    0.0001    0.0108   
                       Call Option Values  0.0985    0.3383    0.5781   
                       Call Option Values  1.0830    1.3306    1.5749   

  Strike=25.00, Sigma= 0.40, Alpha= 0.00
                       Call Option Values  0.0000    0.0000    0.0000   
                       Call Option Values  0.0000   -0.0000    0.0000   
                       Call Option Values  0.0000   -0.0000    0.0000

Example 3

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

\[\begin{split}p(x) = \begin{cases} 0, & x \leq K \\ B, & x > K \\ \end{cases}\end{split}\]

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

\[p(x) = \max\left(x - K_1\right) - \max\left(x - K_2\right), \phantom{...} K_2 > K_1\]

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:

1. Strike and bet prices \(K_1 = {10.0}\), \(K_2 = {15.0}\), and \(B = {2.0}\).
2. Volatility \(\sigma = \{0.4\}\).
3. Times until expiration = \(\{1/4, 1/2\}\).
4. Interest rate \(r = 0.1\).
5. \(x_{min} = 0\), \(x_{max} = 30\).
6. \(nx = 61\), \(n = 3 \times nx = 183\).

from __future__ import print_function
from numpy import *
from pyimsl.math.feynmanKac import feynmanKac
from pyimsl.math.feynmanKacEvaluate import feynmanKacEvaluate

# These routines define the coefficients, payoff, boundary conditions
# and forcing term for American and European Options.


def fkcfiv_call(x, tx, iflag, value):
    if iflag[0] == 0:
        # The payoff function - Use flag passed to decide which
        if type == 1:
            # After reaching the strike price the payoff jumps
            # from zero to the bet value.
            if x > strike_price:
                value[0] = bet
            else:
                value[0] = 0.0
        else:  # type=2
            # Function is zero up to strike price.
            # Then linear between strike price and spread.
            # Then has constant value Spread-Strike Price after
            # the value Spread.
            value[0] = max(x - strike_price, 0.0) - max(x - spread, 0.0)
    elif iflag[0] == -1:
        # The coefficient derivative d sigma/ dx
        value[0] = sigma
    elif iflag[0] == 1:
        # The coefficient sigma(x)
        value[0] = sigma * x
    elif iflag[0] == 2:
        # The coefficient mu(x)
        value[0] = (interest_rate - dividend) * x
    elif iflag[0] == 3:
        # The coefficient kappa(x)
        value[0] = interest_rate

    # Note that there is no time dependence
    iflag[0] = 0


def fkbcp_call(nbc, tx, iflag, val):
    if iflag[0] == 1:
        val[0] = 1.0
        val[1] = 0.0
        val[2] = 0.0
        val[3] = 0.0
        val[4] = 0.0
        val[5] = 1.0
        val[6] = 0.0
        val[7] = 0.0
        val[8] = 0.0
        val[9] = 0.0
        val[10] = 1.0
        val[11] = 0.0
        # Note no time dependence in case (1) for IFLAG
        iflag[0] = 0

    if iflag[0] == 2:
        # This is the discount factor using the risk-free
        # interest rate
        df = exp(interest_rate * tx)
        # Use flag passed to decide on boundary condition
        if type == 1:
            val[0] = 1.0
            val[1] = 0.0
            val[2] = 0.0
            val[3] = bet * df
        else:  # type=2
            val[0] = 1.0
            val[1] = 0.0
            val[2] = 0.0
            val[3] = (spread - strike_price) * df
        val[4] = 0.0
        val[5] = 1.0
        val[6] = 0.0
        val[7] = 0.0
        val[8] = 0.0
        val[9] = 0.0
        val[10] = 1.0
        val[11] = 0.0


NXGRID = 61
NTGRID = 2
NV = 12

# The strike price
strike_price = 10.0

# The spread value
spread = 15.0

# The Bet for the Cash-or-Nothing Call
bet = 2.0

# The sigma value
sigma = 0.4

# Time values for the options
nt = 2
time = [0.25, 0.5]

# Values of the underlying where evaluations are made
xs = empty(NV)

# Value of the interest rate and continuous dividend
interest_rate = 0.1
dividend = 0.0

# Values of the min and max underlying values modeled
x_min = 0.0
x_max = 30.0

# Define parameters for the integration step.
nint = NXGRID - 1
n = 3 * NXGRID

xgrid = empty(NXGRID)
yb = empty(0)
ybprime = empty(0)
yv = empty(0)
yvprime = empty(0)
fb = empty([NTGRID, NV])
fv = empty([NTGRID, NV])

# Number of left/right boundary conditions.
nlbcd = 3
nrbcd = 3

# Define an equally-spaced grid of points for the underlying
# price
dx = (x_max - x_min) / ((float)(nint))
xgrid[0] = x_min
xgrid[NXGRID - 1] = x_max

for i in range(2, NXGRID):
    xgrid[i - 1] = xgrid[i - 2] + dx

for i in range(1, NV + 1):
    xs[i - 1] = 2.0 + (i - 1) * 2.0

# Flag the difference in payoff functions
# 1 for the Bet, 2 for the Vertical Spread
type = 1
feynmanKac(nlbcd, nrbcd, xgrid, time,
           fkcfiv_call, fkbcp_call, yb, ybprime)

type = 2
feynmanKac(nlbcd, nrbcd, xgrid, time,
           fkcfiv_call, fkbcp_call, yv, yvprime)

# Evaluate solutions at vector of points XS(:), at each time value
# prior to expiration.

for i in range(0, NTGRID):
    tmp = feynmanKacEvaluate(xgrid, xs, yb[i, ])
    fb[i, :] = tmp
    tmp = feynmanKacEvaluate(xgrid, xs, yv[i, ])
    fv[i, :] = tmp

print("  European Option Value for A Bet")
print("   and a Vertical Spread, 3 and 6 Months Prior to Expiry")
print("     Number of equally spaced spline knots:%4d" % NXGRID)
print("     Number of unknowns:%4d" % n)
print("     Strike=%5.2f, Sigma=%5.2f, Interest Rate=%5.2f" %
      (strike_price, sigma, interest_rate))
print("     Bet=%5.2f, Spread Value=%5.2f\n" % (bet, spread))
print("       Underlying             A Bet   Vertical Spread")
for i in range(0, NV):
    print("        %9.4f%9.4f%9.4f%9.4f%9.4f" %
          (xs[i], fb[0, i], fb[1, i], fv[0, i], fv[1, i]))

Output

  European Option Value for A Bet
   and a Vertical Spread, 3 and 6 Months Prior to Expiry
     Number of equally spaced spline knots:  61
     Number of unknowns: 183
     Strike=10.00, Sigma= 0.40, Interest Rate= 0.10
     Bet= 2.00, Spread Value=15.00

       Underlying             A Bet   Vertical Spread
           2.0000   0.0000  -0.0000  -0.0000  -0.0000
           4.0000   0.0000   0.0000  -0.0000   0.0000
           6.0000   0.0000   0.0112  -0.0000   0.0038
           8.0000   0.0017   0.2686  -0.0000   0.1486
          10.0000   1.0000   0.9948   0.0124   0.8898
          12.0000   1.9983   1.6103   2.0000   2.1904
          14.0000   2.0000   1.8650   4.0003   3.4267
          16.0000   2.0000   1.9335   5.0003   4.2274
          18.0000   2.0000   1.9477   5.0000   4.6261
          20.0000   2.0000   1.9501   5.0000   4.7903
          22.0000   2.0000   1.9505   5.0000   4.8493
          24.0000   2.0000   1.9506   5.0000   4.8685

Example 4

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 \(\nu x, \nu \geq 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:

1. Bond face value \(K = \{ 1 \}\), conversion factor \(\nu = 1.125\)
2. Volatility \(\sigma = \{ 0.25 \}\)
3. Times until expiration \(= \{ 1/2, 1 \}\)
4. Interest rate \(r = 0.1\) , dividend \(D = 0.02\)
5. \(x_{\min} = 0, \phantom{...} x_{\max} = 4\)
6. \(nx = 61, \phantom{...} n = 3 \times nx = 183\)
7. Boundary conditions \(f(0,t) = K \exp(-r(T-t)), f \left(x_{\max}, t \right) = vx_{\max}\)
8. Terminal data \(f(x,T) = \max(K, \nu x)\)
9. Constraint for bond holder \(f(x,t) \geq \nu x\)

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

from __future__ import print_function
from numpy import *
from pyimsl.math.feynmanKac import feynmanKac
from pyimsl.math.feynmanKacEvaluate import feynmanKacEvaluate

# These routines define the coefficients, payoff, boundary conditions
# and forcing term.


def fkcfiv_cbond(x, tx, iflag, value):
    if iflag[0] == 0:
        # The payoff function
        value[0] = max(factor * x, strike_price)
    elif iflag[0] == -1:
        # The coefficient derivative d sigma/ dx
        value[0] = sigma
    elif iflag[0] == 1:
        # The coefficient sigma(x)
        value[0] = sigma * x
    elif iflag[0] == 2:
        # The coefficient mu(x)
        value[0] = (interest_rate - dividend) * x
    elif iflag[0] == 3:
        # The coefficient kappa(x)
        value[0] = interest_rate

    # Note that there is no time dependence
    iflag[0] = 0


def fkbcp_cbond(nbc, tx, iflag, val):
    if iflag[0] == 1:
        dp = strike_price * exp(tx * interest_rate)
        val[0] = 1.0
        val[1] = 0.0
        val[2] = 0.0
        val[3] = dp
        val[4] = 0.0
        val[5] = 1.0
        val[6] = 0.0
        val[7] = 0.0
        val[8] = 0.0
        val[9] = 0.0
        val[10] = 1.0
        val[11] = 0.0
    elif iflag[0] == 2:
        val[0] = 1.0
        val[1] = 0.0
        val[2] = 0.0
        val[3] = factor * x_max
        val[4] = 0.0
        val[5] = 1.0
        val[6] = 0.0
        val[7] = factor
        val[8] = 0.0
        val[9] = 0.0
        val[10] = 1.0
        val[11] = 0.0
        # Note no time dependence
        iflag[0] = 0


def fkforce_cbond(interval, ndeg, nxgrid, y,
                  time, width, xlocal, qw, u, phi, dphi):
    local = 6
    yl = empty(6)
    bf = empty(6)
    for i in range(0, local):
        yl[i] = y[3 * interval - 3 + i]
        phi[i] = 0.0
    for i in range(0, local * local):
        dphi[i] = 0.0
    mu = 2

    # This is the local definition of the forcing term -
    # It "forces" the constraint f >= factor*x.
    for j in range(1, local + 1):
        for l in range(1, ndeg + 1):
            bf[0] = u[(l - 1)]
            bf[1] = u[(l - 1) + ndeg]
            bf[2] = u[(l - 1) + 2 * ndeg]
            bf[3] = u[(l - 1) + 6 * ndeg]
            bf[4] = u[(l - 1) + 7 * ndeg]
            bf[5] = u[(l - 1) + 8 * ndeg]

            rt = 0.0
            for k in range(0, local):
                rt += yl[k] * bf[k]

            rt = value / (rt + value - factor * xlocal[l - 1])
            phi[j - 1] += qw[l - 1] * bf[j - 1] * pow(rt, mu)

    for i in range(0, local):
        phi[i] = -phi[i] * width * factor * strike_price

    # This is the local derivative matrix for the forcing term
    for j in range(1, local + 1):
        for i in range(1, local + 1):
            for l in range(1, ndeg + 1):
                bf[0] = u[(l - 1)]
                bf[1] = u[(l - 1) + ndeg]
                bf[2] = u[(l - 1) + 2 * ndeg]
                bf[3] = u[(l - 1) + 6 * ndeg]
                bf[4] = u[(l - 1) + 7 * ndeg]
                bf[5] = u[(l - 1) + 8 * ndeg]

                rt = 0.0
                for k in range(0, local):
                    rt += yl[k] * bf[k]

                rt = 1.0 / (rt + value - factor * xlocal[l - 1])
                dphi[i - 1 + (j - 1) * local] += qw[l - 1] * bf[i - 1] * \
                    bf[j - 1] * pow(value * rt, mu) * rt

    for i in range(0, local * local):
        dphi[i] = -mu * dphi[i] * width * factor * strike_price


def fkinit_cbond(nxgrid, ntgrid, xgrid, tgrid,
                 time, yprime, y, atol, rtol):
    if time == 0.0:
        # Set initial data precisely.
        for i in range(1, nxgrid):
            if xgrid[i - 1] * factor < strike_price:
                y[3 * i - 3] = strike_price
                y[3 * i - 2] = 0.0
                y[3 * i - 1] = 0.0
            else:
                y[3 * i - 3] = xgrid[i - 1] * factor
                y[3 * i - 2] = factor
                y[3 * i - 1] = 0.0

# Compute value of a Convertible Bond


NXGRID = 61
NTGRID = 2
NV = 13

# The face value
strike_price = 1.0

# The sigma  or volatility value
sigma = 0.25

# Time values for the options
time = [0.5, 1.0]

# Values of the underlying where evaluation are made
xs = empty(NV)

# Value of the interest rate, continuous dividend and factor
interest_rate = 0.1
dividend = 0.02
factor = 1.125

# Values of the min and max underlying values modeled
x_min = 0.0
x_max = 4.0

# Define parameters for the integration step.
nint = NXGRID - 1
n = 3 * NXGRID

xgrid = empty(NXGRID)
y = empty(0)
yprime = empty(0)
f = empty([NTGRID + 1, NV])

# Number of left/right boundary conditions.
nlbcd = 3
nrbcd = 3

# Define an equally-spaced grid of points for the
# underlying price
dx = (x_max - x_min) / float(nint)
xgrid[0] = x_min
xgrid[NXGRID - 1] = x_max

for i in range(2, NXGRID):
    xgrid[i - 1] = xgrid[i - 2] + dx

for i in range(1, NV + 1):
    xs[i - 1] = (i - 1) * 0.25

# Use a pure absolute error tolerance for the integration
value = 1.0e-3

# Compute value of convertible bond
feynmanKac(nlbcd, nrbcd, xgrid, time,
           fkcfiv_cbond, fkbcp_cbond, y, yprime,
           fcnInit=fkinit_cbond,
           fcnForce=fkforce_cbond,
           atolRtolScalars={'atol': 1.0e-3, 'rtol': 0.0e0})

# Evaluate and display solutions at vector of points XS(:),
# at each time value prior to expiration.
for i in range(0, NTGRID + 1):
    tmp = feynmanKacEvaluate(xgrid, xs, y[i, ])
    f[i, :] = tmp

print("  Convertible Bond Value, 0+, 6 and 12 Months Prior to Expiry")
print("     Number of equally spaced spline knots:%4d" % NXGRID)
print("     Number of unknowns:%4d" % n)
print("     Strike=%5.2f, Sigma=%5.2f" % (strike_price, sigma))
print("     Interest Rate=%5.2f, Dividend=%5.2f, Factor=%6.3f\n" %
      (interest_rate, dividend, factor))
print("     Underlying        Bond Value")

for i in range(0, NV):
    print("       %8.4f%8.4f%8.4f%8.4f" %
          (xs[i], f[0, i], f[1, i], f[2, i]))

Output

  Convertible Bond Value, 0+, 6 and 12 Months Prior to Expiry
     Number of equally spaced spline knots:  61
     Number of unknowns: 183
     Strike= 1.00, Sigma= 0.25
     Interest Rate= 0.10, Dividend= 0.02, Factor= 1.125

     Underlying        Bond Value
         0.0000  1.0000  0.9512  0.9048
         0.2500  1.0000  0.9512  0.9049
         0.5000  1.0000  0.9513  0.9065
         0.7500  1.0000  0.9737  0.9605
         1.0000  1.1250  1.1416  1.1464
         1.2500  1.4062  1.4117  1.4121
         1.5000  1.6875  1.6922  1.6922
         1.7500  1.9688  1.9731  1.9731
         2.0000  2.2500  2.2540  2.2540
         2.2500  2.5312  2.5349  2.5349
         2.5000  2.8125  2.8160  2.8160
         2.7500  3.0938  3.0970  3.0970
         3.0000  3.3750  3.3781  3.3781

Fatal Errors

IMSL_STOP_USER_FCN

Request from user supplied function to stop algorithm. User flag = “#”.