monteCarlo

monteCarlo(evolve, payoff, initialState, timeGrid, numPaths, [numFactors=1], [params=NULL], [options=NULL])

Performs generic Monte Carlo simulation. It generates simulated paths according to a user-defined stochastic process, evaluates a user-defined objective function along the paths, and estimates the expected value of that objective function: θ = E[f(X)]. Where:

• X is a user-defined random variable or stochastic process, defined in this function by evolve;
• f(X) is a user-defined objective function, defined in this function by payoff.

This function generates random drivers, handles correlations among random factors, simulates paths in batches, and aggregates statistical results. It supports single-state or multi-state processes, single-factor or multi-factor drivers, and single-output or multi-output functions. Pricing financial derivatives is one of its primary use cases; users should define product rules and discounting logic in payoff.

evolve is a function definition, which specifies the state evolution of a user-defined stochastic process. The function signature is fixed as: evolve(t, state, dt, noise, params).

Where:

• batchPathCount is the actual number of paths in the current batch;
• stateDim is the state dimension, equal to the number of elements in initialState; a scalar corresponds to 1;
• numFactors is the numFactors parameter passed when calling monteCarlo.

evolve is called once per batch and per time step and returns a DOUBLE matrix of shape batchPathCount × stateDim.

payoff is a function definition, which specifies a user-defined objective function.

• When the "pathMode" of the options parameter is set to "full", the function signature is fixed as: payoff(paths, timeGrid, params).

• When the "pathMode" of the options parameter is set to "terminal", the function signature is fixed as: payoff(terminalState, timeGrid, params). Where:

  • terminalState is a DOUBLE matrix of size batchPathCount × stateDim containing only the terminal states of the current batch;
  • The timeGrid and params parameters are the same as those in payoff(paths, timeGrid, params).

The objective function is called per batch:

• For single output, returns a numeric vector of length batchPathCount, where each element corresponds to an objective function sample for one simulated path;
• For multiple outputs, returns a numeric matrix of shape batchPathCount × numOutputs, where each column corresponds to an independent target to be estimated.

initialState is a numeric scalar or vector that specifies the initial state of the stochastic process, that is, the initial state at the first time point in timeGrid. If specified as a scalar, it represents a single-state process; if specified as a vector, it represents a multi-state process.

timeGrid is a numeric vector that specifies the simulation time grid. Must contain at least two finite non-NULL values and be strictly increasing. Supports non-equally spaced time grids.

numPaths is an INT or LONG scalar that specifies the total number of simulated paths. Must be greater than or equal to 2; when antithetic variates are enabled, it must be an even number greater than 2.

numFactors (optional) is an INT scalar. Specifies the number of random factors input to the state evolution function at each time step. Defaults to 1 and must be greater than 0. The number of random factors is independent of the state dimension.

params (optional) is an object of ANY type. Defaults to NULL. Specifies a user-defined parameter object passed to the evolve and payoff functions. Within a single simulation, callbacks should treat it as a read-only object.

options (optional) is a dictionary (STRING, ANY). Defaults to NULL. Supports the following keys:

The format of the return value is determined by "resultMode" in options:

• When set to "value", returns only the Monte Carlo estimate. A single-output objective function corresponds to a DOUBLE scalar, and a multi-output objective function corresponds to a DOUBLE vector.
• When set to "detail", returns a dictionary of statistical results containing the following keys:

For a single output, "estimate", "variance", "stdError", "confidenceLower", and "confidenceUpper" are scalars; for multiple outputs, they are vectors whose length equals the number of output columns.

Example 1. Basic Pattern of General Pricing

This example demonstrates the basic pattern of general pricing: evolve defines the model, payoff defines the product, and options controls the numerical computation. The results are, in order, the estimate, the standard error, and the lower and upper bounds of the 95% confidence interval.

def gbmEvolve(t, state, dt, noise, p) {
    mu = (p["r"] - p["q"] - 0.5 * p["sigma"] * p["sigma"]) * dt
    return state * exp(mu + p["sigma"] * sqrt(dt) * noise)
}
// Define an appropriate payoff function for your product
def callPayoff(terminalState, timeGrid, p) {
    terminal = terminalState[, 0]
    T = timeGrid[size(timeGrid) - 1] - timeGrid[0]
    return exp(-p["r"] * T) * max(terminal - p["K"], 0.0)
}
p = dict(STRING, DOUBLE)
p["r"] = 0.03
p["q"] = 0.0
p["sigma"] = 0.20
p["K"] = 100.0
mcOpt = dict(STRING, ANY)
mcOpt["seed"] = 42
mcOpt["antithetic"] = true
mcOpt["batchSize"] = 2000
mcOpt["pathMode"] = "terminal"
mcResult = monteCarlo(gbmEvolve, callPayoff, 100.0, 0..12 / 12.0, 20000, 1, p, mcOpt)
[mcResult["estimate"], mcResult["stdError"], mcResult["confidenceLower"], mcResult["confidenceUpper"]]
// output: [9.184365113062904, 0.073151245422059, 9.040991306611417, 9.32773891951439]

Example 2. Valuing a European Call Option Using the Full Path

In this example, you define geometric Brownian motion and the discounted option payoff. Fix the random seed and batch size, and enable antithetic variates. The results are, in order, the estimate, the standard error, and the lower and upper bounds of the 95% confidence interval.

def gbmEvolve(t, state, dt, noise, params) {
    drift = (params["rate"] - 0.5 * params["volatility"] * params["volatility"]) * dt
    return state * exp(drift + params["volatility"] * sqrt(dt) * noise)
}
def callPayoff(paths, timeGrid, params) {
    terminal = paths[size(paths) - 1][, 0]
    maturity = timeGrid[size(timeGrid) - 1] - timeGrid[0]
    return exp(-params["rate"] * maturity) * max(terminal - params["strike"], 0.0)
}
params = dict(keyType=STRING, valueType=DOUBLE)
params["rate"] = 0.03
params["volatility"] = 0.2
params["strike"] = 100.0
options = dict(keyType=STRING, valueType=ANY)
options["seed"] = 42
options["batchSize"] = 1000
options["antithetic"] = true
options["pathMode"] = "full"
result = monteCarlo(evolve=gbmEvolve, payoff=callPayoff, initialState=100.0,
    timeGrid=0..12 / 12.0, numPaths=10000, params=params, options=options)
[result["estimate"], result["stdError"], result["confidenceLower"], result["confidenceUpper"]]
// output: [9.32114319573, 0.106148772868, 9.11309542391, 9.52919096756]

Example 3. Simultaneously Valuing Call, Put, and Digital Options Using Terminal Mode

The objective function returns three sample columns. All three columns come from the same set of paths, so the estimates can be compared directly. Specifying pathMode="terminal" means the engine does not need to store intermediate time points. The results are, in order, three estimates, three standard errors, three analytical values, and the column-wise differences between the two.

def gbmEvolve(t, state, dt, noise, p) {
    drift = (p["r"] - p["q"] - 0.5 * p["sigma"] * p["sigma"]) * dt
    return state * exp(drift + p["sigma"] * sqrt(dt) * noise)
}
def vanillaPayoff(terminalState, timeGrid, p) {
    ST = terminalState[, 0]
    T = timeGrid[size(timeGrid) - 1] - timeGrid[0]
    df = exp(-p["r"] * T)
    callPV = df * max(ST - p["K"], 0.0)
    putPV = df * max(p["K"] - ST, 0.0)
    digitalPV = df * iif(ST > p["K"], 1.0, 0.0)
    return matrix(callPV, putPV, digitalPV)
}
p = dict(STRING, DOUBLE)
p["S0"] = 100.0
p["K"] = 100.0
p["r"] = 0.03
p["q"] = 0.0
p["sigma"] = 0.20
opt = dict(STRING, ANY)
opt["seed"] = 42
opt["antithetic"] = true
opt["batchSize"] = 10000
opt["pathMode"] = "terminal"
result = monteCarlo(gbmEvolve, vanillaPayoff, p["S0"], 0..12 / 12.0, 100000, 1, p, opt)
T = 1.0
d1 = (log(p["S0"] / p["K"]) + (p["r"] - p["q"] + 0.5 * p["sigma"] * p["sigma"]) * T) / (p["sigma"] * sqrt(T))
d2 = d1 - p["sigma"] * sqrt(T)
bsCall = p["S0"] * exp(-p["q"] * T) * cdfNormal(0.0, 1.0, d1) - p["K"] * exp(-p["r"] * T) * cdfNormal(0.0, 1.0, d2)
bsPut = p["K"] * exp(-p["r"] * T) * cdfNormal(0.0, 1.0, -d2) - p["S0"] * exp(-p["q"] * T) * cdfNormal(0.0, 1.0, -d1)
bsDigital = exp(-p["r"] * T) * cdfNormal(0.0, 1.0, d2)
analytic = [bsCall, bsPut, bsDigital]
assert 1, all(abs(result["estimate"] - analytic) <= 4.0 * result["stdError"])
[result["estimate"], result["stdError"], analytic, result["estimate"] - analytic]
// output: [[9.414867268425256, 6.464866000587869, 0.503573891813656], [0.033089109740657, 0.021276101501186, 0.000413951476158], [9.413403383853016, 6.457956738703836, 0.504572291844394], [0.001463884572241, 0.006909261884034, -0.000998400030737]]

Example 4. Asian Call Option (Arithmetic Average and Geometric Average)

The payoff of an Asian option depends on the observed prices along the entire path. This example samples only at the 12 month-end observation points, excludes t=0 from the average, and uses the same set of paths to produce both arithmetic-average and geometric-average estimates. The results are, in order, two estimates and two standard errors.

def gbmEvolve(t, state, dt, noise, p) {
    drift = (p["r"] - p["q"] - 0.5 * p["sigma"] * p["sigma"]) * dt
    return state * exp(drift + p["sigma"] * sqrt(dt) * noise)
}
def asianPayoff(paths, timeGrid, p) {
    n = rows(paths[0])
    m = size(paths) - 1
    // Use only the 12 future observation times, excluding t=0
    prices = matrix(DOUBLE, n, m)
    for (i in 1..m) {
        prices[, i - 1] = paths[i][, 0]
    }
    arithmeticAverage = rowAvg(prices)
    geometricAverage = exp(rowAvg(log(prices)))
    T = timeGrid[size(timeGrid) - 1] - timeGrid[0]
    df = exp(-p["r"] * T)
    return df * matrix(max(arithmeticAverage - p["K"], 0.0), max(geometricAverage - p["K"], 0.0))
}
p = dict(STRING, DOUBLE)
p["r"] = 0.03
p["q"] = 0.0
p["sigma"] = 0.20
p["K"] = 100.0
opt = dict(STRING, ANY)
opt["seed"] = 42
opt["antithetic"] = true
opt["batchSize"] = 5000
opt["pathMode"] = "full"
res = monteCarlo(gbmEvolve, asianPayoff, 100.0, 0..12 / 12.0, 100000, 1, p, opt)
[res["estimate"], res["stdError"]]
// output: [[5.609002269061869, 5.413219715290193], [0.019006421439442, 0.018509710160461]]

Example 5. Seasoned Asian Call Option

The contract has entered its averaging period. At the valuation time, six fixings have already been determined and the remaining six have not yet occurred. Historical accumulated values cannot be written into initialState; instead, they must be passed in with params. The remaining observation times should be located in timeGrid and used from there. The results are, in order, the estimate, the standard error, and the lower and upper bounds of the 95% confidence interval.

def gbmEvolve(t, state, dt, noise, p) {
    drift = (p["r"] - p["q"] - 0.5 * p["sigma"] * p["sigma"]) * dt
    return state * exp(drift + p["sigma"] * sqrt(dt) * noise)
}
def seasonedAsianPayoff(paths, timeGrid, p) {
    n = rows(paths[0])
    futureSum = take(0.0, n)
    futureCount = size(paths) - 1
    // Start from the first pending observation date after the valuation date; do not include that date's price again
    for (i in 1..futureCount) {
        futureSum += paths[i][, 0]
    }
    averagePrice = (p["realizedSum"] + futureSum) / (p["realizedCount"] + futureCount)
    remainingT = timeGrid[size(timeGrid) - 1] - timeGrid[0]
    return exp(-p["r"] * remainingT) * max(averagePrice - p["K"], 0.0)
}
p = dict(STRING, DOUBLE)
p["r"] = 0.03
p["q"] = 0.0
p["sigma"] = 0.20
p["K"] = 100.0
p["realizedCount"] = 6.0
p["realizedSum"] = 98.0 + 99.0 + 101.0 + 102.0 + 96.0 + 104.0
opt = dict(STRING, ANY)
opt["seed"] = 2026
opt["antithetic"] = true
opt["batchSize"] = 5000
opt["pathMode"] = "full"
seasoned = monteCarlo(gbmEvolve, seasonedAsianPayoff, 102.0, 0..6 / 12.0, 100000, 1, p, opt)
[seasoned["estimate"], seasoned["stdError"], seasoned["confidenceLower"], seasoned["confidenceUpper"]]
// output: [2.633867593334059, 0.006737018412437, 2.620663279882499, 2.647071906785618]

Example 6. Two-Asset Basket Call Option

The basket option holds two stocks with equal weights. The returns of the two stocks are correlated; "corr" applies to the random factors rather than to the state variables themselves. The example also demonstrates numFactors=2 and a vector-valued initialState.

def twoAssetEvolve(t, state, dt, noise, p) {
    s1 = state[, 0] * exp((p["rate"] - 0.5 * p["vol1"] * p["vol1"]) * dt + p["vol1"] * sqrt(dt) * noise[, 0])
    s2 = state[, 1] * exp((p["rate"] - 0.5 * p["vol2"] * p["vol2"]) * dt + p["vol2"] * sqrt(dt) * noise[, 1])
    return matrix(s1, s2)
}
def basketPayoff(terminalState, timeGrid, p) {
    basketPrice = 0.5 * terminalState[, 0] + 0.5 * terminalState[, 1]
    maturity = timeGrid[size(timeGrid) - 1] - timeGrid[0]
    return exp(-p["rate"] * maturity) * max(basketPrice - p["strike"], 0.0)
}
p = dict(STRING, ANY)
p["rate"] = 0.03
p["vol1"] = 0.20
p["vol2"] = 0.25
p["strike"] = 100.0
options = dict(STRING, ANY)
options["seed"] = 42
options["antithetic"] = true
options["batchSize"] = 5000
options["pathMode"] = "terminal"
options["corr"] = matrix(1.0 0.4, 0.4 1.0)
result = monteCarlo(twoAssetEvolve, basketPayoff, 100.0 95.0, 0..12 / 12.0, 100000, 2, p, options)
[result["estimate"], result["stdError"]]
// output: [7.52196546515141, 0.030194782803288]

Example 7. Barrier Options and Lookback Options

The payoffs of the up-and-out call and the lookback call both depend on the maximum price along the entire path, so you must use the default pathMode="full". Both share the same set of paths, and a single simulation produces two estimates. Here, whether the barrier is hit is determined only at 12 discrete observation points, so the result cannot be compared directly with the theoretical price under continuous monitoring. The results are, in order, two estimates and two standard errors.

def gbmEvolve(t, state, dt, noise, p) {
    drift = (p["r"] - p["q"] - 0.5 * p["sigma"] * p["sigma"]) * dt
    return state * exp(drift + p["sigma"] * sqrt(dt) * noise)
}
def barrierLookbackPayoff(paths, timeGrid, p) {
    maxSpot = paths[0][, 0]
    for (i in 1..(size(paths) - 1)) {
        maxSpot = max(maxSpot, paths[i][, 0])
    }
    ST = paths[size(paths) - 1][, 0]
    T = timeGrid[size(timeGrid) - 1] - timeGrid[0]
    df = exp(-p["r"] * T)
    hit = maxSpot >= p["barrier"]
    upOut = df * iif(hit, 0.0, max(ST - p["K"], 0.0))
    lookback = df * max(maxSpot - p["K"], 0.0)
    return matrix(upOut, lookback)
}
p = dict(STRING, DOUBLE)
p["r"] = 0.03
p["q"] = 0.0
p["sigma"] = 0.20
p["K"] = 100.0
p["barrier"] = 120.0
opt = dict(STRING, ANY)
opt["seed"] = 42
opt["antithetic"] = true
opt["batchSize"] = 5000
opt["pathMode"] = "full"
res = monteCarlo(gbmEvolve, barrierLookbackPayoff, 100.0, 0..12 / 12.0, 100000, 1, p, opt)
[res["estimate"], res["stdError"]]
// output: [[1.80604717983654, 14.783698314987175], [0.011698546708195, 0.028289631948896]]

Example 8. Two-Asset Basket Option and Exchange Option

The basket option holds two stocks with equal weights. The exchange option exchanges asset 2 for asset 1 at maturity. The prices of the two underlying assets are correlated. This example specifies numFactors=2, a vector-valued initialState, and "corr": the correlation matrix acts only on the random factors and does not directly change the state variables. The results are, in order, two estimates, two standard errors, the Margrabe analytic value, and the difference.

// corr only correlates the random factors; it does not directly modify the state variables
def twoAssetEvolve(t, state, dt, noise, p) {
    s1 = state[, 0] * exp((p["r"] - p["q1"] - 0.5 * p["sig1"] * p["sig1"]) * dt + p["sig1"] * sqrt(dt) * noise[, 0])
    s2 = state[, 1] * exp((p["r"] - p["q2"] - 0.5 * p["sig2"] * p["sig2"]) * dt + p["sig2"] * sqrt(dt) * noise[, 1])
    return matrix(s1, s2)
}
def twoAssetPayoff(terminal, timeGrid, p) {
    s1 = terminal[, 0]
    s2 = terminal[, 1]
    T = timeGrid[size(timeGrid) - 1] - timeGrid[0]
    df = exp(-p["r"] * T)
    basket = df * max(0.5 * s1 + 0.5 * s2 - p["basketK"], 0.0)
    exchange = df * max(s1 - s2, 0.0)
    return matrix(basket, exchange)
}
p = dict(STRING, DOUBLE)
p["S1"] = 100.0
p["S2"] = 95.0
p["r"] = 0.03
p["q1"] = 0.01
p["q2"] = 0.02
p["sig1"] = 0.20
p["sig2"] = 0.25
p["rho"] = 0.40
p["basketK"] = 100.0
R = matrix(1.0 0.0, 0.0 1.0)
R[0, 1] = p["rho"]
R[1, 0] = p["rho"]
opt = dict(STRING, ANY)
opt["seed"] = 42
opt["antithetic"] = true
opt["batchSize"] = 10000
opt["pathMode"] = "terminal"
opt["corr"] = R
res = monteCarlo(twoAssetEvolve, twoAssetPayoff, [p["S1"], p["S2"]], 0..12 / 12.0, 200000, 2, p, opt)
T = 1.0
sigM = sqrt(p["sig1"] * p["sig1"] + p["sig2"] * p["sig2"] - 2.0 * p["rho"] * p["sig1"] * p["sig2"])
d1 = (log(p["S1"] / p["S2"]) + (p["q2"] - p["q1"] + 0.5 * sigM * sigM) * T) / (sigM * sqrt(T))
d2 = d1 - sigM * sqrt(T)
reference = p["S1"] * exp(-p["q1"] * T) * cdfNormal(0.0, 1.0, d1) - p["S2"] * exp(-p["q2"] * T) * cdfNormal(0.0, 1.0, d2)
assert 1, abs(res["estimate"][1] - reference) <= 4.0 * res["stdError"][1]
[res["estimate"], res["stdError"], reference, res["estimate"][1] - reference]
// output: [[6.775993366037865, 12.799129569741291], [0.020845620528879, 0.020807928733168], 12.783129860200672, 0.015999709540617]

Example 9. Heston Stochastic Volatility Model

Variance itself is a stochastic process. The state is two-dimensional, (S, v), and requires two random factors with a correlation of -0.7: noise[, 0] drives the price, and noise[, 1] drives the variance. The variance uses a full-truncation Euler update and is truncated to zero both before and after the update to avoid negative variance. The second output column is the expected variance at maturity, not an option cash flow, so it is not discounted. The results are, in order, two estimates and two standard errors.

// state has two columns: [stock price S, instantaneous variance v]
// noise has two columns: two correlated standard normal drivers
def hestonEvolve(t, state, dt, noise, p) {
    S = state[, 0]
    v = state[, 1]
    vPos = max(v, 0.0)
    nextS = S * exp((p["r"] - p["q"] - 0.5 * vPos) * dt + sqrt(vPos * dt) * noise[, 0])
    nextV = max(v + p["kappa"] * (p["theta"] - vPos) * dt + p["xi"] * sqrt(vPos * dt) * noise[, 1], 0.0)
    return matrix(nextS, nextV)
}
def hestonPayoff(terminalState, timeGrid, p) {
    T = timeGrid[size(timeGrid) - 1] - timeGrid[0]
    ST = terminalState[, 0]
    vT = terminalState[, 1]
    callPV = exp(-p["r"] * T) * max(ST - p["K"], 0.0)
    // The second column is E[vT], not an option cash flow, and must not be discounted again
    return matrix(callPV, vT)
}
p = dict(STRING, DOUBLE)
p["r"] = 0.03
p["q"] = 0.0
p["K"] = 100.0
p["kappa"] = 2.0
p["theta"] = 0.04
p["xi"] = 0.50
p["rho"] = -0.70
R = matrix(1.0 0.0, 0.0 1.0)
R[0, 1] = p["rho"]
R[1, 0] = p["rho"]
opt = dict(STRING, ANY)
opt["seed"] = 42
opt["batchSize"] = 10000
opt["antithetic"] = true
opt["pathMode"] = "terminal"
opt["corr"] = R
res = monteCarlo(hestonEvolve, hestonPayoff, 100.0 0.04, 0..252 / 252.0, 300000, 2, p, opt)
[res["estimate"], res["stdError"]]
// output: [[8.945060830037322, 0.040530651223046], [0.012914704996569, 6.5769063872e-05]]

Example 10. Snowball Autocallable Structure

The payoff at maturity depends on whether the knock-out price is hit and whether the knock-in price is hit during the life of the contract; you cannot look only at the terminal value. After the first knock-out, that path generates no further cash flows, and the payment must be discounted over the actual remaining life rather than always to maturity. This example checks for knock-out every 63 simulation days (approximately quarterly) and records knock-in on every simulation day. It specifies pathMode="full" to advance the path state point by point. The payment is based on a notional principal of 100, so the result is on the order of 100 and is not an option price. The results are, in order, the present value of the payment, the standard error, and the lower and upper bounds of the 95% confidence interval.

def gbmEvolve(t, state, dt, noise, p) {
    drift = (p["r"] - p["q"] - 0.5 * p["sigma"] * p["sigma"]) * dt
    return state * exp(drift + p["sigma"] * sqrt(dt) * noise)
}
def snowballPayoff(paths, timeGrid, p) {
    n = rows(paths[0])
    alive = take(true, n)
    knockedIn = take(false, n)
    pv = take(0.0, n)
    last = size(paths) - 1
    for (i in 1..last) {
        spot = paths[i][, 0]
        t = timeGrid[i] - timeGrid[0]
        // Check for knock-in on every simulation day; the state after knock-out does not affect subsequent payments
        knockedIn = knockedIn || (alive && (spot <= p["kiLevel"]))
        // In this example, simulation days 63, 126, 189, and 252 are knock-out observation days
        if (i == 63 || i == 126 || i == 189 || i == 252) {
            newKO = alive && (spot >= p["koLevel"])
            redemption = p["notional"] * (1.0 + p["coupon"] * t)
            pv += iif(newKO, exp(-p["r"] * t) * redemption, 0.0)
            alive = alive && (!newKO)
        }
    }
    ST = paths[last][, 0]
    T = timeGrid[last] - timeGrid[0]
    noKnockInPay = p["notional"] * (1.0 + p["coupon"] * T)
    knockInPay = p["notional"] * (1.0 + min(ST / p["referenceSpot"] - 1.0, 0.0))
    maturityPay = iif(knockedIn, knockInPay, noKnockInPay)
    // Only paths that are not knocked out early receive the cash flow at maturity
    pv += iif(alive, exp(-p["r"] * T) * maturityPay, 0.0)
    return pv
}
p = dict(STRING, DOUBLE)
p["r"] = 0.03
p["q"] = 0.0
p["sigma"] = 0.20
p["notional"] = 100.0
p["referenceSpot"] = 100.0
p["koLevel"] = 103.0
p["kiLevel"] = 80.0
p["coupon"] = 0.12
opt = dict(STRING, ANY)
opt["seed"] = 42
opt["antithetic"] = true
opt["batchSize"] = 5000
opt["pathMode"] = "full"
res = monteCarlo(gbmEvolve, snowballPayoff, 100.0, 0..252 / 252.0, 100000, 1, p, opt)
[res["estimate"], res["stdError"], res["confidenceLower"], res["confidenceUpper"]]
// output: [99.09039224470335, 0.036261836916808, 99.01932035033315, 99.16146413907356]

Example 11. Digital Options and Range Binary Options

The payoff is a discontinuous 0/1 jump, expressed using iif. Digital options are far more sensitive to the number of paths than European options, and typically require more samples to achieve the same accuracy.

def gbmEvolve(t, state, dt, noise, p) {
    return state * exp((p["rate"] - 0.5 * p["volatility"] * p["volatility"]) * dt + p["volatility"] * sqrt(dt) * noise)
}
def digitalPayoff(terminalState, timeGrid, p) {
    maturity = timeGrid[size(timeGrid) - 1] - timeGrid[0]
    discount = exp(-p["rate"] * maturity)
    spot = terminalState[, 0]
    cashDigital = discount * iif(spot > p["strike"], 1.0, 0.0)
    rangeBinary = discount * iif((spot >= p["lower"]) && (spot <= p["upper"]), 1.0, 0.0)
    return matrix(cashDigital, rangeBinary)
}
p = dict(STRING, ANY)
p["rate"] = 0.03
p["volatility"] = 0.2
p["strike"] = 100.0
p["lower"] = 90.0
p["upper"] = 110.0
options = dict(STRING, ANY)
options["seed"] = 42
options["antithetic"] = true
options["batchSize"] = 5000
options["pathMode"] = "terminal"
result = monteCarlo(gbmEvolve, digitalPayoff, 100.0, 0..12 / 12.0, 200000, 1, p, options)
[result["estimate"], result["stdError"]]
// output: [[0.504947072243628, 0.372025147513488], [0.000303013624505, 0.001406220595986]]

Example 12. Retrieving Full Paths for Distribution Analysis

Beyond pricing, you may also want to examine the terminal price distribution (for risk reports, stress tests, or plotting). Set returnPaths=true to include full paths in the returned result, so you can compute user-defined statistics such as quantiles directly from the paths. Note that it is available only with the default pathMode="full" and resultMode="detail".

def gbmEvolve(t, state, dt, noise, p) {
    return state * exp((p["rate"] - 0.5 * p["volatility"] * p["volatility"]) * dt + p["volatility"] * sqrt(dt) * noise)
}
def terminalPayoff(paths, timeGrid, p) {
    return paths[size(paths) - 1][, 0]
}
p = dict(STRING, ANY)
p["rate"] = 0.03
p["volatility"] = 0.2
options = dict(STRING, ANY)
options["seed"] = 42
options["batchSize"] = 5000
options["returnPaths"] = true
result = monteCarlo(gbmEvolve, terminalPayoff, 100.0, 0..12 / 12.0, 20000, 1, p, options)
terminalPrices = result["paths"][size(result["paths"]) - 1][, 0]
[size(result["paths"]), avg(terminalPrices), percentile(terminalPrices, 5.0), percentile(terminalPrices, 95.0)]
// output: [13, 102.8214156238338, 72.57171400837616, 139.8270247863367]

Example 13. Interest Rate Caplet

The payoff of an interest rate product requires both the terminal interest rate and the discount factor. Accumulate the discount factor as a second state along the interest rate path so that it can be multiplied directly at the valuation time, avoiding repeated integration in payoff.

def vasicekEvolve(t, state, dt, noise, p) {
    rate = state[, 0]
    nextRate = rate + p["kappa"] * (p["theta"] - rate) * dt + p["sigma"] * sqrt(dt) * noise[, 0]
    nextDiscount = state[, 1] * exp(-0.5 * (rate + nextRate) * dt)
    return matrix(nextRate, nextDiscount)
}
def capletPayoff(terminalState, timeGrid, p) {
    terminalRate = terminalState[, 0]
    terminalDiscount = terminalState[, 1]
    return terminalDiscount * max(terminalRate - p["strike"], 0.0) * p["accrual"]
}
p = dict(STRING, ANY)
p["kappa"] = 0.8
p["theta"] = 0.03
p["sigma"] = 0.01
p["strike"] = 0.025
p["accrual"] = 0.25
options = dict(STRING, ANY)
options["seed"] = 42
options["antithetic"] = true
options["batchSize"] = 5000
options["pathMode"] = "terminal"
result = monteCarlo(vasicekEvolve, capletPayoff, 0.03 1.0, 0..12 / 12.0, 200000, 1, p, options)
[result["estimate"], result["stdError"]]
// output: [0.001462868116922, 1.236656302e-06]

Example 14. Antithetic Variates

When you set antithetic=true, Z and -Z must be paired, and both numPaths and batchSize must be even; otherwise, the call is rejected before execution:

def gbmEvolve(t, state, dt, noise, p) {
    return state * exp((p["rate"] - 0.5 * p["volatility"] * p["volatility"]) * dt + p["volatility"] * sqrt(dt) * noise)
}
def callPayoff(terminalState, timeGrid, p) {
    maturity = timeGrid[size(timeGrid) - 1] - timeGrid[0]
    return exp(-p["rate"] * maturity) * max(terminalState[, 0] - p["strike"], 0.0)
}
p = dict(STRING, ANY)
p["rate"] = 0.03
p["volatility"] = 0.2
p["strike"] = 100.0
options = dict(STRING, ANY)
options["seed"] = 42
options["antithetic"] = true
options["batchSize"] = 5000
options["pathMode"] = "terminal"
monteCarlo(gbmEvolve, callPayoff, 100.0, 0..12 / 12.0, 10001, 1, p, options)
// output: antithetic simulation requires even numPaths and batchSize.

A correct example is shown below:

def gbmEvolve(t, state, dt, noise, p) {
    return state * exp((p["rate"] - 0.5 * p["volatility"] * p["volatility"]) * dt + p["volatility"] * sqrt(dt) * noise)
}
def callPayoff(terminalState, timeGrid, p) {
    maturity = timeGrid[size(timeGrid) - 1] - timeGrid[0]
    return exp(-p["rate"] * maturity) * max(terminalState[, 0] - p["strike"], 0.0)
}
p = dict(STRING, ANY)
p["rate"] = 0.03
p["volatility"] = 0.2
p["strike"] = 100.0
options = dict(STRING, ANY)
options["seed"] = 42
options["antithetic"] = true
options["batchSize"] = 5000
options["pathMode"] = "terminal"
r = monteCarlo(gbmEvolve, callPayoff, 100.0, 0..12 / 12.0, 10000, 1, p, options)
plainOptions = dict(STRING, ANY)
plainOptions["seed"] = 42
plainOptions["antithetic"] = false
plainOptions["batchSize"] = 5000
plainOptions["pathMode"] = "terminal"
plain = monteCarlo(gbmEvolve, callPayoff, 100.0, 0..12 / 12.0, 10000, 1, p, plainOptions)
[r["estimate"], r["stdError"], plain["estimate"], plain["stdError"], plain["stdError"] / r["stdError"]]
// output: [9.350067041131123, 0.105390522690999, 9.279319321832183, 0.138518059304004, 1.314331267813654]

Related Functions: mcSimulate, mcEstimate