monteCarlo

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

执行通用蒙特卡洛数值求解，按用户定义的随机过程生成模拟路径，在路径上计算用户定义的目标函数，并给出该目标函数期望值的估计 θ = E[f(X)]。其中：

• X 为用户定义的随机变量或随机过程，本函数中由 evolve 定义；

• f(X) 为用户定义的目标函数，本函数中由 payoff 定义。

本函数负责生成随机驱动、处理随机因子相关性、按批次模拟路径并聚合统计结果，支持单状态或多状态过程、单因子或多因子驱动，以及单输出或多输出目标函数。金融衍生品定价是其主要应用场景之一，产品规则和折现逻辑由用户在 payoff 中定义。

evolve 函数定义，指定自定义随机过程状态演化，函数签名固定为：evolve(t, state, dt, noise, params)。

其中：

• batchPathCount 为当前批次实际包含的路径数；

• stateDim 为状态维度，等于 initialState 的元素个数，标量对应 1；

• numFactors 为调用 monteCarlo 传入的 numFactors 参数。

evolve 每批、每个时间步调用一次，返回 batchPathCount × stateDim 的 DOUBLE 矩阵。

payoff 函数定义，指定自定义目标函数。

• 当指定 options 参数的 "pathMode" 为 "full" 时，函数签名固定为：payoff(paths, timeGrid, params)。

• 当指定 options 参数的 "pathMode" 为 "terminal" 时，函数签名固定为：payoff(terminalState, timeGrid, params)。其中：
  • terminalState 是 batchPathCount × stateDim 的 DOUBLE 矩阵，仅包含当前批次的终值状态；
  • timeGrid 与 params 参数与 payoff(paths, timeGrid, params) 一致。

目标函数按批次调用：

• 单输出时，返回长度为 batchPathCount 的数值向量，每个元素对应一条模拟路径的目标函数样本；
• 多输出时，返回形为 batchPathCount × numOutputs 的数值矩阵，每列对应一个独立待估计目标。

initialState 数值标量或向量，指定随机过程的初始状态，即 timeGrid 第一个时间点的初始状态。指定为标量时表示单状态过程，指定为向量则表示多状态过程。

timeGrid 数值向量，指定模拟时间网格。至少包含两个非 NULL 的有限数值，且严格递增。支持非等间隔时间网格。

numPaths INT 或 LONG 标量，指定模拟路径总数。必须大于等于 2；启用对偶变量时，必须为大于 2 的偶数。

numFactors 可选参数， INT 标量，指定每个时间步输入状态演化函数的随机因子数，默认值为 1，必须大于 0。随机因子数与状态维度相互独立。

params 可选参数， ANY 类型对象，默认值为 NULL。指定用户自定义参数对象，传递给 evolve 和 payoff 函数。同一次模拟中的回调应将其作为只读对象使用。

options 可选参数，字典（STRING, ANY），默认值为 NULL。支持以下键：

返回值形式由 options 中的 "resultMode" 决定：

• 指定为 "value" 时，仅返回 Monte Carlo 估计值。单输出目标函数对应 DOUBLE 标量，多输出目标函数对应 DOUBLE 向量。
• 指定为 "detail" 时，返回统计结果字典，包含以下键：

单输出时，"estimate"、"variance"、"stdError"、"confidenceLower"、"confidenceUpper" 为标量；多输出时，为与输出列数相同的向量。

例 1. 通用定价的基本模式

本例演示通用定价的基本模式：evolve 定义模型，payoff 定义产品，options 控制数值计算。结果依次为估值、标准误、95% 置信区间下限和上限。

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)
}
// 按自己的产品定义合适的收益函数
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]

例 2. 使用完整路径计算欧式看涨期权的估值

本例由用户定义几何布朗运动及折现后的期权回报。固定随机种子和批次大小，开启对偶变量。结果依次为估值、标准误、95% 置信区间下限和上限。

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]

例 3. 使用终值模式同时计算看涨、看跌与数字期权的估值

目标函数返回三列样本，三列来自同一批路径，估计值之间可直接比较；指定 pathMode="terminal" 使引擎不必保存中间时点。结果依次为三个估计值、三个标准误、三个解析值，以及两者的逐列差值。

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]]

例 4. 亚式看涨期权（算术平均与几何平均）

亚式期权的赔付依赖整条路径上的观察价格。本例只在 12 个月末观察点取样、不把 t=0 计入平均值，并用同一批路径同时给出算术平均与几何平均两个估计值。结果依次为两个估计值和两个标准误。

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
    // 只用未来 12 个观察时点，不含 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]]

例 5. 存续期亚式看涨期权

合约已进入平均期，估值时点前已有 6 个 Fixing 已确定，剩余 6 个尚未发生。历史累计值不能写入 initialState，而应随 params 传入；剩余观察时点在 timeGrid 中定位后取用。结果依次为估值、标准误、95% 置信区间下限和上限。

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
    // 从估值日后的第一个待观察日开始，不重复计入当日价格
    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]

例 6. 双资产篮子看涨期权

篮子期权按等权持有两只股票，两者收益相关，"corr" 作用于随机因子而非状态本身。示例同时演示了 numFactors=2 与向量型 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]

例 7. 障碍期权与回望期权

向上敲出看涨与回望看涨的赔付都依赖整条路径上的最高价，必须用默认的 pathMode="full"。两者共用同一批路径，一次模拟给出两个估计值。这里只按 12 个离散观察点判断障碍是否被触及，与连续监测的理论价格不可直接比较。结果依次为两个估计值和两个标准误。

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]]

例 8. 双资产篮子期权与交换期权

篮子期权按等权持有两只股票，交换期权到期以资产 2 换资产 1，两个标的价格相关。本例指定 numFactors=2、向量型 initialState 与 "corr"：相关矩阵只作用于随机因子，不直接改变状态变量。结果依次为两个估计值、两个标准误、Margrabe 解析值和差值。

// corr 仅相关化随机因子，不直接修改状态变量
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]

例 9. Heston 随机波动率模型

方差本身是随机过程，状态为 (S, v) 两维，并需要两个相关系数为 -0.7 的随机因子：noise[, 0] 驱动价格，noise[, 1] 驱动方差。方差采用全截断 Euler 更新，更新前后都截断到零，避免出现负方差。第二列输出是到期方差的期望，不是期权现金流，因此不参与折现。结果依次为两个估计值和两个标准误。

// state 两列：[股价 S, 瞬时方差 v]
// noise 两列：已相关化的两个标准正态驱动
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)
    // 第二列是 E[vT]，不是期权现金流，不能再折现
    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]]

例 10. 雪球自动赎回结构

到期赔付取决于存续期内是否触及敲出价、是否触及敲入价，不能只看终值。首次敲出后该路径不再产生未来现金流，且兑付须按实际存续时间折现，不能一律按到期日折现。本例按每 63 个模拟日（约每季度）观察敲出、每个模拟日记录敲入，指定 pathMode="full" 逐点推进路径状态。兑付按名义本金 100 计，因此结果量级约为 100，不是期权价格。结果依次为兑付现值、标准误、95% 置信区间下限和上限。

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]
        // 每个模拟日观察敲入，敲出后的状态不影响后续给付
        knockedIn = knockedIn || (alive && (spot <= p["kiLevel"]))
        // 本例第 63、126、189、252 个模拟日为敲出观察日
        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)
    // 仅未提前敲出的路径领取到期现金流
    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]

例 11. 数字期权与区间二元期权

赔付是不连续的 0/1 跳变，靠 iif 表达。数字期权对路径数远比欧式期权敏感，同一精度下通常需要更多样本。

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]]

例 12. 取回完整路径做分布分析

定价之外还想看期末价格分布（风险报告、压力测试、画图）。指定 returnPaths=true 使返回结果额外带完整路径，可以直接对路径做分位数等自定义统计。注意它仅在默认的 pathMode="full" 且 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]

例 13. 利率上限单元

利率产品的赔付要同时用到期末利率和贴现因子。把贴现因子作为第二个状态随利率路径累积，回到估值时点直接相乘，避免在 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]

例 14. 对偶变量法

指定 antithetic=true 时， Z 与 -Z 必须成对，numPaths 与 batchSize 都须为偶数，否则会在执行前被拒绝：

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.

正确示例如下：

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]

相关函数： mcSimulate, mcEstimate