Finite Difference Method for MultiAsset Option Valuation
Description
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.
Usage
1 
Arguments
X 
list of inputs. List items given in the Details section. 
Details
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.
Value
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.
Author(s)
Michael Eichenberger and Carlo Rosa
References
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.
Examples
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)
