FEYNMAN_KAC Function
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 C2 Hermite quintic splines.
Usage
result = FEYNMAN_KAC (nlbcd, nrbcd, xgrid, tgrid, fcn_fkcfiv, fcn_fkbcp, y, y_prime)
Input Parameters
nlbcd—A scalar long value indicating the number of left boundary conditions. It is required that 1 ≤ nlbcd ≤ 3.
nrbcd—A scalar long value indicating the number of right boundary conditions. It is required that 1 ≤ nrbcd ≤ 3.
xgrid—Array of length nxgrid containing the breakpoints for the Hermite quintic splines used in the x discretization, where nxgrid is the number of grid lines in the x-direction. The points in xgrid must be strictly increasing. The values xgrid(0) and xgrid(nxgrid – 1) are the endpoints of the interval.
tgrid—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.
fcn_fkcfiv—User-supplied function, fcn_fkcfiv ( x, t, iflag, value), to compute the values of the coefficients σ, σ', μ, κ for the Feynman-Kac PDE and the initial data function p(x), xmin ≤ x ≤ xmax.
x—Space variable.
t—Time variable.
iflag—On entry, iflag indicates which coefficient or data function is to be computed. The following table shows which value has to be returned by fcn_fkcfiv for all possible values of iflag:
iflag | Computed coefficient | |
|---|---|---|
–1 |  | |
0 | p(x) | |
1 | σ | |
2 | μ | |
3 | κ |
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.
value—Value of the coefficient or initial data function. Which value is computed depends on the input value for iflag.
fcn_fkbcp—User-supplied function, fcn_fkbcp ( nbc, t, iflag, values), to define boundary values that the solution of the differential equation must satisfy. There are nlbcd conditions specified at the left end, xmin, and nrbcd conditions at the right end, xmax. The boundary conditions are:
nbc—Number of boundary conditions. At xmin, nbc=nlbcd, at xmax, nbc = nrbcd.
t—Time point of the boundary conditions.
iflag—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 = xmin | |
2 | Right end, x = xmax |
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.
values—Array of length 4*max(nlbcd, nrbcd) containing the values of the boundary condition coefficients in its first 4*nbc locations. The coefficients for xmin are stored row-wise according to the following scheme:
The coefficients for xmax are stored similarly.
Output Parameters
y—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:
for t = tgrid(i–1), i = 1, ..., ntgrid. The representation for the initial data at t = 0 is:
The (ntgrid+1) by (3*nxgrid) matrix:
is stored row-wise in array y.
After the integration, use row i of array y as input argument coef in function FEYNMAN_KAC_EVALUATE Function to obtain an array of values for f(x, t) or its partials fx, fxx, fxxx at time point t=0, tgrid(i–1), i=1, ..., ntgrid.
The expressions for the basis functions βj(x) are represented piece-wise and can be found in Hanson, R. (2008) “Integrating Feynman-Kac Equations Using Hermite Quintic Finite Elements”.
y_prime—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.
for 
and:
for 
The (ntgrid + 1) by (3*nxgrid) matrix:
is stored row-wise in array y_prime.
After the integration, use row i of array y_prime as input argument coef in function FEYNMAN_KAC_EVALUATE to obtain an array of values for the partials ft, ftx, ftxx, ftxxx 
at time point t=tgrid(i–1), i=1, ..., ntgrid and row 0 for the partials at t=0.
at time point t=tgrid(i–1), i=1, ..., ntgrid and row 0 for the partials at t=0.Returned Value
result—This function returns 1 to indicate success.
Input Keywords
Double—If present and nonzero, then double precision is used.
Fcn_fkinit—User-supplied function, Fcn_fkinit(nxgrid, ntgrid, xgrid, tgrid, time, yprime, y, atol, rtol), 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, fcn_fkinit can change the contents of y to match data for any reason.
nxgrid—Number of grid lines in the x-direction.
ntgrid—Number of time points where an approximate solution is returned.
xgrid—Vector of length nxgrid containing the breakpoints for the Hermite quintic splines used in the x discretization.
tgrid—Vector of length ntgrid containing the set of time points (in time-remaining units) where an approximate solution is returned.
time—Time variable.
yprime—Vector of length 3*nxgrid containing the derivative of the coefficients of the Hermite quintic spline at time point time.
y—Vector of length 3*nxgrid containing the coefficients of the Hermite quintic spline at time point time.
atol—Array of length 3*nxgrid containing absolute error tolerances.
rtol—Array of length 3*nxgrid containing relative error tolerances.
Fcn_force—User-supplied function, Fcn_force(interval, ndeg, nxgrid, y, time, width, xlocal, qw, u, phi, dphi), used in case there is a non-zero term φ(f,x,t) in the Feynman-Kac differential equation. If function Fcn_force is not used, it is assumed that φ(f,x,t) is identically zero.
interval—Index indicating the bounds xgrid(interval–1) and xgrid(interval) of the integration interval, 1 ≤ interval ≤ nxgrid–1.
ndeg—The degree used for the Gauss-Legendre formulas.
nxgrid—Number of grid lines in the x-direction.
y—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.
time—Time variable.
width—The interval width, width = xgrid(interval) – xgrid(interval–1).
xlocal—Array of length ndeg containing the Gauss-Legendre points translated and normalized to the interval [xgrid(interval–1), xgrid(interval)].
qw—Vector of length ndeg containing the Gauss-Legendre weights.
u—Array of dimension 12 by ndeg containing the basis function values that define 
at the Gauss-Legendre points xlocal. Setting
at the Gauss-Legendre points xlocal. Setting
and 
, 
is defined as:
phi—Vector of length 6 containing Gauss-Legendre approximations for the local contributions:
where t = time and:
Vector phi contains elements:
for i=0, ..., 5.
dphi—Array of dimension 6 by 6 containing a Gauss-Legendre approximation for the Jacobian of the local contributions ϕt at t = time:
The approximation to this symmetric matrix is stored row-wise, i.e.:
for i,j=0, ..., 5.
ScalarAtolRtol—Scalar values that apply to the error estimates of all components of the solution y in the differential equation solver SDASLX. See the Atol and Rtol keywords 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.
Atol—Componentwise tolerances are used for the computation of solution y in the differential equation solver SDASLX. Atol is an array 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 keyword 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.
Rtol—Componentwise tolerances are used for the computation of solution y in the differential equation solver SDASLX. Rtol is an array 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 keyword 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.
Ndeg—The degree used for the Gauss-Legendre formulas for constructing the finite element matrices. It is required that Ndeg ≥ 6. Default: Ndeg = 6.
Tdepend—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 | time dependence of variable |
|---|---|
0 | σ' |
1 | σ |
2 | μ |
3 | κ |
4 | Left boundary conditions |
5 | Right boundary conditions |
6 | φ |
Max_stepsize—This is the maximum step size the integrator may take. Default: Max_stepsize = MACHINE(MAX_POS), the largest possible machine number.
Init_stepsize—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: Init_stepsize = –ε, where ε is the machine precision
Max_Nsteps—The maximum number of steps between each output point of the integration. Default: Max_Nsteps = 500000.
Step_control—Indicates which step control algorithm is used for the integration. If Step_control = 0, then the step control method of Söderlind is used. If Step_control = 1, then the method used by the original Petzold code SASSL is used. Default: Step_control = 0.
Max_bdf_order—Maximum order of the backward differentiation formulas used in the integrator. It is required that 1 ≤ Max_bdf_order ≤ 5. Default: Max_bdf_order = 5.
T_barrier—This keyword controls whether the code should integrate past a special point, T_barrier, and then interpolate to get y and y' at the points in tgrid. If this keyword is present, the integrator assumes the equations either change on the alternate sides of T_barrier or they are undefined there. In this case, the code creeps up to T_barrier in the direction of integration. It is required that T_barrier ≥ tgrid(ntgrid–1). Default: T_barrier = tgrid(ntgrid–1).
Output Keywords
Istate—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 condition data |
5 | Right boundary condition data |
6 | Forcing term |
Default: Istate(i) = 0 for i = 0, ..., 6.
Nval—Array of length 3 summarizing the number of evaluations required during the integration.
i | nval(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. |
Discussion
The generalized Feynman-Kac differential equation has the form
where the initial data satisfies f(x, T) = p(x). The derivatives are 
, etc. The FEYNMAN_KAC function uses a finite element Galerkin method over the rectangle:
, etc. The FEYNMAN_KAC function uses a finite element Galerkin method over the rectangle: 
in (x, t) to compute the approximate solution. The interval [xmin, xmax] is decomposed with a grid:
(xmin =)x1 < x2 < ... < xm (= xmax)
On each subinterval the solution is represented by:
The values:
are time-dependent coefficients associated with each interval. The basis functions b0, b1, b2 are given for:
by:
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:
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 FEYNMAN_KAC_EVALUATE Function.
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 σ = {0.4}
3. Times until expiration = {1/4, 1/2}
4. Interest rate r = 0.1
5. xmin = 0.0, xmax = 30.0
6. nx = 61, n = 3 × nx = 183
The payoff function is the “vanilla option”, p(x) = max(K – x, 0).
; These routines define the coefficients, payoff, boundary
; conditions and forcing term for American and European
; Options.
FUNCTION fkcfiv_put, x, tx, iflag, value
; The sigma value.
sigma=0.4
; Value of the interest rate and continuous dividend.
strike_price=10.0
interest_rate=0.1
dividend=0.0
zero=0.0
CASE iflag OF
0:BEGIN
; The payoff function.
value = MAX([strike_price - x, zero])
END
-1:BEGIN
; The coefficient derivative dsigma/dx.
value = sigma
END
1:BEGIN
; The coefficient sigma(x).
value = sigma*x
END
2:BEGIN
; The coefficient mu(x).
value = (interest_rate - dividend) * x
END
3:BEGIN
; The coefficient kappa(x).
value = interest_rate
END
ENDCASE
; Note that there is no time dependence.
iflag = 0L
RETURN,1
END
FUNCTION fkbcp_put, nbc, tx, iflag, val
CASE iflag OF
1:BEGIN
val(*) = 0.0
val(1) = 1.0
val(3) = -1.0
val(6) = 1.0
END
2:BEGIN
val(*) = 0.0
val(0) = 1.0
val(5) = 1.0
val(10) = 1.0
END
ENDCASE
; Note no time dependence.
iflag = 0
RETURN,1
END
FUNCTION fkforce_put, interval, ndeg, nxgrid, $
y, time, width, xlocal, $
qw, u, phi, dphi
COMMON f_kac_usr, usr_data
local=6
yl = FLTARR(local)
bf = FLTARR(local)
zero=0.0
one=1.0
index = INDGEN(local)
data = usr_data
yl = y(3*interval-3 + index)
phi(*) = zero
strike_price = data(0)
interest_rate = data(1)
value = data(2)
mu=2.0
; This is the local definition of the forcing term.
FOR j=1L, local DO BEGIN
FOR l=1L, ndeg DO BEGIN
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 = SUM(yl * bf)
rt = value/(rt + value - (strike_price-xlocal(l-1)))
phi(j-1) = phi(j-1) + (qw(l-1) * bf(j-1) * rt^mu)
ENDFOR
ENDFOR
phi = -phi * width*interest_rate*strike_price
; This is the local derivative matrix for the forcing term.
dphi(*) = zero
FOR j=1L, local DO BEGIN
FOR i=1L, local DO BEGIN
FOR l=1L, ndeg DO BEGIN
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 = SUM(yl * bf)
rt = one/(rt + value-(strike_price-xlocal(l-1)))
dphi(i-1+(j-1)*local) = dphi(i-1+(j-1)*local) + $
qw(l-1) * bf(i-1) * bf(j-1) $
* (rt^(mu+1.0))
ENDFOR
ENDFOR
ENDFOR
FOR i=0L, (local*local)-1 DO BEGIN
dphi(i) = mu * dphi(i) * width * value^mu * $
interest_rate * strike_price
ENDFOR
RETURN,1
END
PRO feynman_kac_ex1
; Compute American Option Premium for Vanilla Put.
COMMON f_kac_usr, usr_data
NXGRD = 61
NTGRD = 2
NV = 9
; The strike price.
KS = 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.
r = 0.1
; Values of the min and max underlying values modeled.
x_min=0.0
x_max=30.0
; Define parameters for the integration step.
nint=NXGRD-1
n=3*NXGRD
xgrid = FLTARR(NXGRD)
ye = FLTARR((NTGRD+1),3*NXGRD)
yeprime = FLTARR((NTGRD+1),3*NXGRD)
ya = FLTARR((NTGRD+1),3*NXGRD)
yaprime = FLTARR((NTGRD+1),3*NXGRD)
fe = FLTARR(NTGRD*NV)
fa = FLTARR(NTGRD*NV)
dx = 0.0
nlbcd = 2
nrbcd = 3
atol = 0.2e-2
; Array for user-defined data.
usr_data = FLTARR(3)
; Define an equally-spaced grid of points for the
; underlying price.
;
; Determine the increment between points (to be used later...)
dx = (x_max-x_min)/ nint
; Create the grid.
xgrid = INTERPOL([x_min,x_max], NXGRD)
Atol_Rtol = [ 0.5e-2, 0.5e-2 ]
; Set up the user data in the COMMON block for the user
; function(s).
usr_data = [KS, r, atol]
ret = FEYNMAN_KAC(nlbcd, nrbcd, xgrid, time, $
"fkcfiv_put", "fkbcp_put", ye, yeprime, $
ScalarAtolRtol=Atol_Rtol)
ret = FEYNMAN_KAC(nlbcd, nrbcd, xgrid, time, $
"fkcfiv_put", "fkbcp_put", ya, $
yaprime, Fcn_force="fkforce_put" , $
ScalarAtolRtol=Atol_Rtol)
; Evaluate solutions at vector of points XS(:), at each
; time value prior to expiration.
FOR i=0L, nt-1 DO BEGIN
fe(i*NV:(i*NV)+NV-1) = FEYNMAN_KAC_EVALUATE(xgrid, $
xs, ye((i+1),*))
fa(i*NV:(i*NV)+NV-1) = FEYNMAN_KAC_EVALUATE(xgrid, xs, $
ya((i+1),*))
ENDFOR
PRINT,'American Option Premium for Vanilla Put, ' + $
'3 and 6 Months Prior to Expiry'
PRINT," Number of equally spaced spline knots: ", NXGRD
PRINT," Number of unknowns: ", n
PRINT," Strike=",KS," sigma=",sigma," Interest Rate=",r
PRINT,""
PRINT,"Underlying European American"
FOR i=0L, NV-1 DO BEGIN
PRINT, xs(i), fe(i), fe(i+NV), fa(i), fa(i+NV), $
format='(5(f10.4))'
ENDFOR
PRINT,""
END
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.0000, sigma= 0.400000, Interest Rate= 0.100000
Underlying European American
0.0000 9.7536 9.5134 10.0000 10.0000
2.0000 7.7537 7.5135 8.0000 8.0000
4.0000 5.7538 5.5150 6.0000 6.0000
6.0000 3.7619 3.5664 4.0000 4.0000
8.0000 1.9073 1.9170 2.0146 2.0880
10.0000 0.6518 0.8579 0.6768 0.9042
12.0000 0.1628 0.3392 0.1679 0.3526
14.0000 0.0374 0.1271 0.0374 0.1321
16.0000 0.0091 0.0477 0.0085 0.0499
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”:
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:
1. power α = {2.0, 1.0, 0.0}
2. strike price K = {15.0, 20.0, 25.0}
3. volatility σ = {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. xmin = 0, xmax = 60
8. nx = 121, n = 3 × nx = 363
With this model the Feynman-Kac differential equation is defined by identifying:
x: S
μ(x,t): rx
κ(x,t): r
φ(f,x,t) ≡ 0
The payoff function is the “vanilla option”, p(x) = max(x – K, 0).
FUNCTION fcn_fkcfiv, x, tx, iflag, value
COMMON f_kac_usr, usr_data
zero=0.0
half=0.5
data = usr_data
strike_price = data(0)
sigma = data(2)
alpha = data(3)
interest_rate = data(4)
dividend = data(5)
CASE iflag OF
0:BEGIN
; The payoff function.
value = (x - strike_price)>zero
END
-1:BEGIN
; The coefficient derivative d sigma/ dx.
value = half * alpha * sigma * x^(alpha*half-1.0)
END
1:BEGIN
; The coefficient sigma(x).
value = sigma * x^(alpha*half)
END
2:BEGIN
; The coefficient mu(x).
value = (interest_rate - dividend) * x
END
3:BEGIN
; The coefficient kappa(x).
value = interest_rate
END
ENDCASE
; Note that there is no time dependence.
iflag = 0
RETURN,1
END
FUNCTION fcn_fkbcp, nbc, tx, iflag, val
COMMON f_kac_usr, usr_data
data = usr_data
strike_price = data(0)
x_max = data(1)
interest_rate = data(4)
CASE iflag OF
1:BEGIN
val(*) = 0.0
val(0) = 1.0
val(5) = 1.0
val(10) = 1.0
; Note no time dependence at left end.
iflag = 0
END
2:BEGIN
df = EXP(interest_rate*tx)
val(*) = 0.0
val(0) = 1.0
val(3) = x_max - df*strike_price
val(5) = 1.0
val(7) = 1.0
val(10) = 1.0
END
ENDCASE
RETURN,1
END
PRO t_feynman_kac_ex2
COMMON f_kac_usr, usr_data
NXGRD = 121l
NTGRD = 3l
NV = 3l
; Compute Constant Elasticity of Variance Model for Vanilla Call
; The set of strike prices.
KS= [ 15.0, 20.0, 25.0]
; The set of sigma values.
sigma = [ 0.2, 0.3, 0.4]
; The set of model diffusion powers.
alpha = [2.0, 1.0, 0.0]
; Time values for the options.
nt = NTGRD
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.
r=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 = NXGRD-1
n = 3*NXGRD
xgrid = FLTARR(NXGRD)
y = FLTARR((NTGRD+1)*3*NXGRD)
yprime = FLTARR((NTGRD+1)*3*NXGRD)
f = FLTARR(NTGRD*NV)
dx = 0.0
; Number of left/right boundary conditions.
nlbcd = 3L
nrbcd = 3L
usr_data = FLTARR(6)
; Define equally-spaced grid of points for the underlying price.
xgrid = INTERPOL([x_min,x_max], NXGRD)
PRINT,"Constant Elasticity of Variance Model for Vanilla Call"
PRINT,"Interest Rate: ",r," Continuous Dividend: ", dividend
PRINT,"Minimum and Maximum Prices of Underlying: ", $
STRING(x_min,Format="(f6.2)"), " ",$
STRING(x_max ,Format="(f6.2)")
PRINT,"Number of equally spaced spline knots:",NXGRD-1
PRINT,"Number of unknowns:",n
PRINT,""
PRINT,"Time in Years Prior to Expiration: ", $
STRING(time(0),Format="(f7.4)"), " ",$
STRING(time(1),Format="(f7.4)"), " ",$
STRING(time(2),Format="(f7.4)")
PRINT,"Option valued at Underlying Prices: ", $
STRING(xs(0),Format="(f5.2)"), " ",$
STRING(xs(1),Format="(f5.2)"), " ",$
STRING(xs(2),Format="(f5.2)")
PRINT,""
FOR i1=1L, 3 DO BEGIN ; Loop over power.
FOR i2=1L, 3 DO BEGIN ; Loop over volatility.
FOR i3=1L, 3 DO BEGIN ; Loop over strike price.
; Pass data through into evaluation routines.
usr_data(0) = KS(i3-1)
usr_data(1) = x_max
usr_data(2) = sigma(i2-1)
usr_data(3) = alpha(i1-1)
usr_data(4) = r
usr_data(5) = dividend
ret = FEYNMAN_KAC(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=0L, NTGRD-1 DO $
f(i*NV:(i*NV)+NV-1) = $
FEYNMAN_KAC_EVALUATE(xgrid, xs, REFORM(y((i+1),*)))
outstr = "Strike="+STRING(KS(i3-1),Format="(f6.2)")+ $
" Sigma="+STRING(sigma(i2-1),Format="(f6.2)")+ $
" Alpha="+STRING(alpha(i1-1),Format="(f6.2)")
PRINT, outstr
FOR i=0L, NV-1 DO BEGIN
outstr = " Call Option Values "+ $
STRING(f(i),Format="(f7.4)") + " "+ $
STRING(f(NV+i),Format="(f7.4)") + " "+ $
STRING(f(2*NV+i),Format="(f7.4)")
PRINT, outstr
ENDFOR
ENDFOR
ENDFOR
ENDFOR
END
Output
Constant Elasticity of Variance Model for Vanilla Call
Interest Rate: 0.0500000 Continuous Dividend: 0.00000
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.2575 4.4729
Call Option Values 5.0624 5.2506 5.4490
Call Option Values 6.0624 6.2486 6.4385
Strike= 20.00 Sigma= 0.20 Alpha= 2.00
Call Option Values 0.1310 0.5956 0.9698
Call Option Values 0.5029 1.0889 1.5100
Call Option Values 1.1980 1.7484 2.1751
Strike= 25.00 Sigma= 0.20 Alpha= 2.00
Call Option Values 0.0000 0.0113 0.0744
Call Option Values 0.0000 0.0373 0.1615
Call Option Values 0.0007 0.1026 0.3132
Strike= 15.00 Sigma= 0.30 Alpha= 2.00
Call Option Values 4.0639 4.3397 4.6618
Call Option Values 5.0626 5.2945 5.5783
Call Option Values 6.0624 6.2709 6.5238
Strike= 20.00 Sigma= 0.30 Alpha= 2.00
Call Option Values 0.3109 1.0275 1.5500
Call Option Values 0.7323 1.5422 2.1024
Call Option Values 1.3765 2.1690 2.7385
Strike= 25.00 Sigma= 0.30 Alpha= 2.00
Call Option Values 0.0006 0.1111 0.3542
Call Option Values 0.0037 0.2169 0.5547
Call Option Values 0.0182 0.3857 0.8222
Strike= 15.00 Sigma= 0.40 Alpha= 2.00
Call Option Values 4.0755 4.5138 4.9674
Call Option Values 5.0661 5.4201 5.8324
Call Option Values 6.0634 6.3578 6.7298
Strike= 20.00 Sigma= 0.40 Alpha= 2.00
Call Option Values 0.5116 1.4645 2.1292
Call Option Values 0.9623 1.9957 2.6951
Call Option Values 1.5815 2.6110 3.3236
Strike= 25.00 Sigma= 0.40 Alpha= 2.00
Call Option Values 0.0083 0.3284 0.7785
Call Option Values 0.0284 0.5165 1.0654
Call Option Values 0.0812 0.7686 1.4104
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.0226 0.1049
Call Option Values 0.1495 0.4107 0.6485
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.0016 0.0785 0.2204
Call Option Values 0.1979 0.4997 0.7541
Call Option Values 1.0836 1.3443 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.0005
Strike= 15.00 Sigma= 0.40 Alpha= 1.00
Call Option Values 4.0624 4.2477 4.4310
Call Option Values 5.0624 5.2477 5.4309
Call Option Values 6.0624 6.2477 6.4309
Strike= 20.00 Sigma= 0.40 Alpha= 1.00
Call Option Values 0.0084 0.1542 0.3443
Call Option Values 0.2482 0.5943 0.8729
Call Option Values 1.0871 1.3790 1.6584
Strike= 25.00 Sigma= 0.40 Alpha= 1.00
Call Option Values 0.0000 0.0000 0.0001
Call Option Values 0.0000 0.0000 0.0008
Call Option Values 0.0000 0.0004 0.0066
Strike= 15.00 Sigma= 0.20 Alpha= 0.00
Call Option Values 4.0627 4.2478 4.4312
Call Option Values 5.0623 5.2479 5.4311
Call Option Values 6.0623 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.0818 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.0625 4.2479 4.4312
Call Option Values 5.0623 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.0029
Call Option Values 0.0895 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.0623 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.0111
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:
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(x – K1) – max(x – K2), K2 > K1
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 K1={10.0}, K2 = {15.0}, and B = {2.0}
2. Volatility σ = {0.4}
3. Times until expiration = {1/4, 1/2}
4. Interest rate r = 0.1.
5. xmin = 0, xmax = 30
6. nx = 61, n = 3 × nx = 183
PRO t_feynman_kac_ex3
COMMON f_kac_usr, usr_data
NXGRD = 61l
NTGRD = 2l
NV = 12l
; The strike price.
KS1 = 10.0
; The spread value.
KS2 = 15.0
; The Bet for the Cash-or-Nothing Call.
bet = 2.0
; The sigma value.
sigma = 0.4
; Time values for the options.
time = [0.25, 0.5]
; Values of the underlying where evaluations are made.
xs = FLTARR(NV)
; Value of the interest rate and continuous dividend.
r = 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 = NXGRD-1
n = 3*NXGRD
xgrid = FLTARR(NXGRD)
yb = FLTARR((NTGRD+1)*3*NXGRD)
ybprime = FLTARR((NTGRD+1)*3*NXGRD)
yv = FLTARR((NTGRD+1)*3*NXGRD)
yvprime = FLTARR((NTGRD+1)*3*NXGRD)
fb = FLTARR(NTGRD*NV)
fv = FLTARR(NTGRD*NV)
; Number of left/right boundary conditions.
nlbcd = 3l
nrbcd = 3l
; Structure for the evaluation routines.
usr_data = { , $idope:I32ARR(1), $
rdope:FLTARR(7) $
}
; Define an equally-spaced grid of points for the underlying
; price.
xgrid = INTERPOL([x_min, x_max], NXGRD)
tmp = FINDGEN(NV)
xs = 2.0+(tmp*2.0)
usr_data.rdope(0) = KS1
usr_data.rdope(1) = bet
usr_data.rdope(2) = KS2
usr_data.rdope(3) = x_max
usr_data.rdope(4) = sigma
usr_data.rdope(5) = r
usr_data.rdope(6) = dividend
; Flag the difference in payoff functions
; 1 for the Bet, 2 for the Vertical Spread.
usr_data.idope(0) = 1
ret = FEYNMAN_KAC(nlbcd, nrbcd, $
xgrid, time, "fkcfiv_call", $
"fkbcp_call", yb, ybprime)
usr_data.idope(0) = 2
ret = FEYNMAN_KAC(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=0L, NTGRD-1 DO BEGIN
fb(i*NV:(i*NV)+NV-1) = FEYNMAN_KAC_EVALUATE(xgrid, $
xs, REFORM(yb((i+1),*)))
fv(i*NV:(i*NV)+NV-1) = FEYNMAN_KAC_EVALUATE(xgrid, $
xs, REFORM(yv((i+1),*)))
ENDFOR
PRINT,""
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: ", $
STRTRIM(NXGRD,2)
PRINT," Number of unknowns: ",STRTRIM(n,2)
PRINT," Strike=", STRING(KS1,Format="(f5.2)"),$
" Sigma=",STRING(sigma,Format="(f5.2)"),$
" Interest Rate=",STRING(r,Format="(f5.2)")
PRINT," Bet=",STRING(bet,Format="(f5.2)"), $
" Spread Value= ",STRING(KS2,Format="(f5.2)")
PRINT,""
PRINT," Underlying A Bet Vertical Spread
FOR i=0L, NV-1 DO BEGIN
PRINT," ",STRING(xs(i),Format="(f9.4)")," ", $
STRING(fb(i),Format="(f9.4)"), " ",$
STRING(fb(i+NV),Format="(f9.4)")," ",$
STRING(fv(i),Format="(f9.4)"), " ",$
STRING(fv(i+NV),Format="(f9.4)")
ENDFOR
END
; These routines define the coefficients, payoff, boundary
; conditions and forcing term for American and European Options.
FUNCTION fkcfiv_call, x, tx, iflag, value
COMMON f_kac_usr, usr_data
zero=0.0
data = usr_data
data_real = data.rdope
data_int = data.idope
strike_price = data_real(0)
bet = data_real(1)
spread = data_real(2)
sigma = data_real(4)
interest_rate = data_real(5)
dividend = data_real(6)
CASE iflag OF
0:BEGIN
; The payoff function - Use flag passed to decide which.
CASE data_int(0) OF
1:BEGIN
; After reaching the strike price the payoff
; jumps from zero to the bet value.
value = zero
IF (x GT strike_price) THEN value = bet
END
2:BEGIN
; 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 = ((x-strike_price)>zero) - (x-spread>zero)
END
ENDCASE
END
-1:BEGIN
; The coefficient derivative d sigma/ dx.
value = sigma
END
1:BEGIN
; The coefficient sigma(x).
value = sigma*x
END
2:BEGIN
; The coefficient mu(x).
value = (interest_rate - dividend)*x
END
3:BEGIN
; The coefficient kappa(x).
value = interest_rate
END
ENDCASE
; Note that there is no time dependence.
iflag = 0
RETURN,1
END
FUNCTION fkbcp_call, nbc, tx, iflag, val
COMMON f_kac_usr, usr_data
data = usr_data
data_int = data.idope
data_real = data.rdope
strike_price = data_real(0)
bet = data_real(1)
spread = data_real(2)
interest_rate = data_real(5)
CASE iflag OF
1:BEGIN
val(*) = 0.0
val(0) = 1.0
val(5) = 1.0
val(10) = 1.0
; Note no time dependence in case (1) for IFLAG.
iflag = 0
END
2:BEGIN
; This is the discount factor using the risk-free
; interest rate.
df = EXP(interest_rate*tx)
; Use flag passed to decide on boundary condition.
CASE data_int(0) OF
1:BEGIN
val(*) = 0.0
val(0) = 1.0
val(3) = bet*df
END
2:BEGIN
val(*) = 0.0
val(0) = 1.0
val(3) = (spread-strike_price)*df
END
ENDCASE
val(5) = 1.0
val(10) = 1.0
END
ENDCASE
RETURN,1
END
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.0014 0.0000 0.0006 6.0000 0.0110 0.0722 0.0039 0.0447 8.0000 0.2691 0.4305 0.1479 0.3831 10.0000 0.9948 0.9781 0.8909 1.1927 12.0000 1.6094 1.4287 2.1911 2.2273 14.0000 1.8655 1.6923 3.4254 3.1552 16.0000 1.9337 1.8177 4.2263 3.8263 18.0000 1.9476 1.8700 4.6264 4.2492 20.0000 1.9501 1.8903 4.7911 4.4921 22.0000 1.9505 1.8979 4.8497 4.6231 24.0000 1.9506 1.9007 4.8684 4.6909Example 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 vx, v ≥ 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 v = 1.125
2. Volatility σ = {0.25}
3. Times until expiration = {1/2, 1}
4. Interest rate r = 0.1, dividend D = 0.02
5. xmin = 0, xmax = 4
6. nx = 61, n = 3 × nx = 183
7. Boundary conditions f(0, t) = K exp(–r(T – t)), f(xmax, t) = vxmax
8. Terminal data f(x, T) = max(K, vx)
9. Constraint for bond holder f(x, t) ≥ vx
Note that f(x, t) ≥ vx the error tolerance is set to a pure absolute error of value 10–3. The free boundary constraint f(x, t) ≥ vx is achieved by use of a non-linear forcing term in the user function fkforce_cbond. The terminal conditions are provided with the user function fkinit_cbond.
PRO t_feynman_kac_ex4
COMMON f_kac_usr, usr_data
; Compute value of a Convertible Bond.
NXGRD = 61
NTGRD = 2
NV = 13
; The face value.
KS = 1.0e0
; The sigma or volatility value.
sigma = 0.25e0
; Time values for the options.
time = [ 0.5, 1.0 ]
; Values of the underlying where evaluation are made.
xs = FLTARR(NV)
; Value of the interest rate, continuous dividend and factor.
r = 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 = NXGRD-1
n = 3*NXGRD
xgrid = FLTARR(NXGRD)
y = FLTARR((NTGRD+1)*3*NXGRD)
yprime = FLTARR((NTGRD+1)*3*NXGRD)
f = FLTARR((NTGRD+1)*NV)
; Array for user-defined data.
usr_data = FLTARR(8)
; Number of left/right boundary conditions.
nlbcd = 3l
nrbcd = 3l
; Define an equally-spaced grid of points for the
; underlying price.
dx = (x_max-x_min)/(nint)
xgrid = INTERPOL([x_min, x_max], NXGRD)
tmp = FINDGEN(NV)
xs = tmp*0.25
; Pass the data for evaluation.
usr_data(0) = KS
usr_data(1) = x_max
usr_data(2) = sigma
usr_data(3) = r
usr_data(4) = dividend
usr_data(5) = factor
; Use a pure absolute error tolerance for the integration.
atol = 1.0e-3
usr_data(6) = atol
; Compute value of convertible bond.
ret = FEYNMAN_KAC(nlbcd, nrbcd, xgrid, time, $
"fkcfiv_cbond", "fkbcp_cbond", y, yprime, $
Fcn_fkinit="fkinit_cbond", $
Fcn_force="fkforce_cbond", $
ScalarAtolRtol=[1.0e-3, 0.0e0] )
; Evaluate and display solutions at vector of points XS(:),
; at each time value prior to expiration.
FOR i=0L, NTGRD DO BEGIN
f(i*NV:(i*NV)+NV-1) = FEYNMAN_KAC_EVALUATE (xgrid, xs, REFORM(y(i,*)))
ENDFOR
PRINT,"Convertible Bond Value, 0+, 6 and 12 Months Prior"+ $
" to Expiry"
PRINT," Number of equally spaced spline knots: ", $
STRTRIM(NXGRD,2)
PRINT," Number of unknowns:",STRTRIM(n,2)
PRINT," Strike= ",STRING(KS,Format="(f5.2)"), $
" Sigma= ",STRING(sigma,Format="(f5.2)")
PRINT," Interest Rate= ",STRING(r,Format="(f5.2)"), $
" Dividend= ",STRING(dividend,Format="(f5.2)"), $
" Factor= ",STRING(factor,Format="(f6.3)")
PRINT,""
PRINT,"Underlying Bond Value"
FOR i=0L, NV-1 DO BEGIN
PRINT," ",STRING(xs(i),Format="(f8.4)")," ", $
STRING(f(i),Format="(f8.4)")," ", $
STRING(f(i+NV),Format="(f8.4)")," ", $
STRING(f(i+2*NV),Format="(f8.4)")
ENDFOR
END
; These routines define the coefficients, payoff,
; boundary conditions and forcing term.
FUNCTION fkcfiv_cbond, x, tx, iflag, value
COMMON f_kac_usr, usr_data
zero=0.0
data = usr_data
strike_price = data(0)
sigma = data(2)
interest_rate = data(3)
dividend = data(4)
factor = data(5)
CASE iflag OF
0:BEGIN
; The payoff function.
value = (factor * x)>strike_price
END
-1:BEGIN
; The coefficient derivative d sigma/ dx.
value = sigma
END
1:BEGIN
; The coefficient sigma(x).
value = sigma*x
END
2:BEGIN
; The coefficient mu(x).
value = (interest_rate - dividend) * x
END
3:BEGIN
; The coefficient kappa(x).
value = interest_rate
END
ENDCASE
; Note that there is no time dependence.
iflag = 0
RETURN,1
END
FUNCTION fkbcp_cbond, nbc, tx, iflag, val
COMMON f_kac_usr, usr_data
data = usr_data
CASE iflag OF
1:BEGIN
strike_price = data(0)
interest_rate = data(3)
dp = strike_price * EXP(tx*interest_rate)
val(*) = 0.0
val(0) = 1.0
val(3) = dp
val(5) = 1.0
val(10) = 1.0
END
2:BEGIN
x_max = data(1)
factor = data(5)
val(*) = 0.0
val(0) = 1.0
val(3) = factor*x_max
val(5) = 1.0
val(7) = factor
val(10) = 1.0
; Note no time dependence.
iflag = 0
END
ENDCASE
RETURN,1
END
FUNCTION fkforce_cbond, interval, ndeg, nxgrid, y, $
time, width, xlocal, qw, $
u, phi, dphi
COMMON f_kac_usr, usr_data
local=6
yl = FLTARR(local)
bf = FLTARR(local)
zero=0.e0
one=1.0e0
data = usr_data
phi(*) = zero
dphi(*) = zero
yl = y((3*interval)-3+INDGEN(local))
value = data(6)
strike_price = data(0)
interest_rate = data(3)
factor = data(5)
mu = 2.0
; This is the local definition of the forcing term.
; It "forces" the constraint f >= factor*x.
FOR j=1L, local DO BEGIN
FOR l=1L, ndeg DO BEGIN
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 = SUM(yl * bf)
rt = value/(rt + value - factor * xlocal(l-1))
phi(j-1) = phi(j-1) + qw(l-1) * bf(j-1) * rt^mu
ENDFOR
ENDFOR
phi = (-phi)*width*factor*strike_price
; This is the local derivative matrix for the forcing term
FOR j=1L, local DO BEGIN
FOR i=1L, local DO BEGIN
FOR l=1L, ndeg DO BEGIN
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 = SUM(yl * bf)
rt = one/(rt + value - factor * xlocal(l-1))
dphi(i-1+(j-1)*local) = $
dphi(i-1+(j-1)*local)+qw(l-1) * bf(i-1) * $
bf(j-1) * ((value*rt)^mu) * rt
ENDFOR
ENDFOR
ENDFOR
dphi = dphi *(-mu) * width * factor * strike_price
RETURN,1
END
FUNCTION fkinit_cbond, nxgrid, ntgrid, xgrid, tgrid, $
time, yprime, y, atol, rtol
COMMON f_kac_usr, usr_data
data = usr_data
IF time EQ 0.0 THEN BEGIN
; Set initial data precisely.
FOR i=1L, nxgrid DO BEGIN
IF (xgrid(i-1) * data(5)) LT data(0) THEN BEGIN
y(3*i-3) = data(0)
y(3*i-2) = 0.0
y(3*i-1) = 0.0
ENDIF ELSE BEGIN
y(3*i-3) = xgrid(i-1) * data(5)
y(3*i-2) = data(5)
y(3*i-1) = 0.0
ENDELSE
ENDFOR
ENDIF
RETURN,1
END
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
