multiAssetOption
generalizes the standard finite difference method to handle mulitple underlying assets, nonuniform grid spacing, nonuniform timestepping, and American exercise. The implementation allows users to vary the option setup (number of underlying assets, call vs. put, European vs. American, etc.) as well as the features of the numerical method (grid spacing, timestepping scheme, etc.). Strike shifting the mesh and Rannacher smoothing are optionally included to remedy problems arising from potential spurious oscillations in the solution.
1 
X 
list of inputs. List items given in the Details section. 
Items of the input list X
are as follows:
X$opt$nAsset
integer; number of underlying assets.
X$opt$payType
case; if 0, digital payoff, if 1, bestof payoff, if 2, worstof payoff.
X$opt$exerType
case; if 0, European exercise, if 1, American exercise.
X$opt$pcFlag
case vector; if 0, call, if 1, put.
X$opt$ttm
scalar; time to maturity.
X$opt$strike
vector; option strikes.
X$opt$rf
scalar; applicable riskfree rate (domestic riskfree rate).
X$opt$q
vector; holding costs of the option's underlyings (dividends, foreign riskfree rates, etc.).
X$opt$vol
vector; volatilities of the option's underlyings.
X$opt$rho
matrix; correlation matrix of the option's underlyings.
X$fd$m
vector; number of spatial steps for each underlying's domain discretization.
X$fd$leftBound
vector; near spatial boundaries of each underlying's domain.
X$fd$kMult
vector; right boundary strike multiples. If 0, far domain boundary calculated via formula given in Kangro and Nicolaides (2000). Otherwise, far domain boundary calculated as the strike multiplied by the strike multiple.
X$fd$density
vector; impacts the concentration of nodes around the option strike. At 0, nodes are uniformly distributed between the near and far boundaries. Increasing the parameter increases the node concentration around the strike.
X$fd$kShift
case vector; if 0, no mesh shifting, if 1, adjusts the node spacing such that the strike falls exactly between two nodes, if 2, adjusts the node spacing such that the strike falls exactly on a node. See Tavella and Randall (2000).
X$fd$theta
scalar; implicitness parameter of the theta method. Chosen between 0 (fully explicit) and 1 (fully implicit).
X$fd$maxSmooth
integer; number of Rannacher smoothing steps. See Rannacher (1984).
X$fd$tol
scalar; error tolerance in penalty iteration for American exercise.
X$fd$maxIter
integer; maximum number of iterations per penalty loop for American exercise.
X$time$tsType
case; if 0, constant timestep size, if 1, adaptive timestep size. See Forsyth and Vetzal (2002).
X$time$N
integer; number of total timesteps if not using adaptive timesteps.
X$time$dtInit
scalar; inital timestep size for adaptive timesteps.
X$time$dNorm
scalar; target relative change for adaptive timesteps.
X$time$D
scalar; normalizing parameter for adaptive timesteps.
The classical order for the state vectors output from the function is illustrated by example. With two underlying assets, option values in each state vector are stored in the order: [11, 21, 31, ... , M1, 12, 22, ... , MN] with M being the total number of nodes used in the first asset spatial discretization and N being the total number of nodes in the second.
multiAssetOption
returns a list:
value 
matrix of perunit option values. Each column stores the state of the option value array (collection of option values for all nodes of the spatial mesh) as a vector following the classical order (see Details section). The columns of the matrix are indexed over time, with the first column beginning at option maturity, and subsequent columns moving backward in time. 
S 
list containing the vectors of spatial grid points associated with each underlying. Vector sizes of underlying spatial grid points need not be equal. 
dimS 
dimension of option value array. This item can be used to reshape the column vectors in 
time 
vector of times associated with each column of the 
For each column (time) of the value
item, the option value at that time can be calculated as the option's notional amount multiplied by the unit option value interpolated over the S
item at the current underlying prices.
Michael Eichenberger and Carlo Rosa
Forsyth, P.A., Vetzal, K.R., 2002. Quadratic convergence for valuing American options using a penalty method. SIAM Journal on Scientific Computing, 23 (6), 2095–2122.
Kangro, R., Nicolaides, R., 2000. Far field boundary conditions for BlackScholes equations. SIAM Journal on Numerical Analysis, 38 (4), 1357–1368.
Rannacher, R., 1984. Finite element solution of diffusion problems with irregular data. Numberische Mathematik, 43, 309–327.
Tavella, D., Randall, C., 2000. Pricing Financial Instruments: The Finite Difference Method. John Wiley & Sons, Inc., New York.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37  # european dualasset digital option example
# initialize inputs list
X < list()
# option inputs
X$opt$nAsset < 2
X$opt$payType < 0
X$opt$exerType < 0
X$opt$pcFlag < c(0, 0)
X$opt$ttm < 0.5
X$opt$strike < c(110, 90)
X$opt$rf < 0.10
X$opt$q < c(0.05, 0.04)
X$opt$vol < c(0.20, 0.25)
X$opt$rho < matrix(c(1, 0.5, 0.5, 1), X$opt$nAsset, X$opt$nAsset)
# finite difference inputs
X$fd$m < c(10, 10)
X$fd$leftBound < c(0, 0)
X$fd$kMult < c(0, 0)
X$fd$density < c(5, 5)
X$fd$kShift < c(1, 1)
X$fd$theta < 0.5
X$fd$maxSmooth < 2
X$fd$tol < 1e7
X$fd$maxIter < 3
# timestep inputs
X$time$tsType < 0
X$time$N < min(X$fd$m) * 4
X$time$dtInit < 0.1 / 4^log2(min(X$fd$m)/5)
X$time$dNorm < 5 / 2^log2(min(X$fd$m)/5)
X$time$D < 0.05
# function check
output < multiAssetOption(X)

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.