mcSimulate

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

Generates full simulation paths using a user-defined state evolution function. The function generates standard normal pseudo-random drivers, handles correlation among random factors, and advances the state along the time grid. It supports multiple states, multiple factors, antithetic variates, and batch execution.

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 to mcSimulate.

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

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 a positive INT or LONG scalar that specifies the total number of simulation paths. 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 that is passed to the evolve function. 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:

Returns a dictionary containing the full simulation paths and simulation information:

Example 1. Generate Geometric Brownian Motion Paths and Random Drivers

This example simulates 8 price paths using 12 equal-length time steps and specifies returnRandomNumbers = true to retrieve the random driver for each step.

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)
}
p = dict(STRING, DOUBLE)
p["r"] = 0.03
p["q"] = 0.0
p["sigma"] = 0.20
opt = dict(STRING, ANY)
opt["seed"] = 7
opt["batchSize"] = 8
opt["returnRandomNumbers"] = true
sim = mcSimulate(gbmEvolve, 100.0, 0..12 / 12.0, 8, 1, p, opt)
[size(sim["paths"]), shape(sim["paths"][0]), size(sim["randomNumbers"]), flatten(sim["paths"][12][, 0])]
// output: [13, [8, 1], 12, [86.07027033006177, 95.76439369271459, 70.83725029862133, 124.81301557000361, 99.78994602443282, 131.41412372603514, 90.4951917379586, 86.17040506214356]]

Example 2. Generate Antithetic Random Drivers for Two Correlated Factors

The correlation matrix specifies a correlation coefficient of 0.5 between factors, using two states and two factors. The first two rows of the returned random drivers form an antithetic pair, and summing them factor by factor yields 0.

def factorEvolve(t, state, dt, noise, params) {
    return state + sqrt(dt) * noise
}
options = dict(keyType=STRING, valueType=ANY)
options["seed"] = 11
options["batchSize"] = 6
options["antithetic"] = true
options["corr"] = matrix(1.0 0.5, 0.5 1.0)
options["returnRandomNumbers"] = true
result = mcSimulate(evolve=factorEvolve, initialState=0.0 0.0,
    timeGrid=0.0 1.0, numPaths=6, numFactors=2, options=options)
noise = result["randomNumbers"][0]
[noise[0, 0] + noise[1, 0], noise[0, 1] + noise[1, 1]]
// output: [0, 0]

Example 3. Quote Multiple Strike Prices in a Single Simulation

The quote sheet requires three strike prices to be quoted at the same time. Calling monteCarlo repeatedly would repeat the simulation, whereas mcSimulate simulates only once, transforms the terminal states into a three-column sample, and passes it to mcEstimate for aggregation—the three quotes are naturally based on the same set of paths, and the covariance between the estimates can be used directly.

def gbmEvolve(t, state, dt, noise, p) {
    return state * exp((p["rate"] - 0.5 * p["volatility"] * p["volatility"]) * dt + p["volatility"] * sqrt(dt) * noise)
}
p = dict(STRING, ANY)
p["rate"] = 0.03
p["volatility"] = 0.2
options = dict(STRING, ANY)
options["seed"] = 42
options["batchSize"] = 5000
simulation = mcSimulate(gbmEvolve, 100.0, 0..12 / 12.0, 50000, 1, p, options)
terminal = simulation["paths"][size(simulation["paths"]) - 1][, 0]
discount = exp(-p["rate"])
samples = discount * matrix(max(terminal - 90.0, 0.0), max(terminal - 100.0, 0.0), max(terminal - 110.0, 0.0))
estimate = mcEstimate(samples)
[estimate["estimate"], estimate["stdError"]]
// output: [[15.411808417905474, 9.397343043691492, 5.253868707329098], [0.075506895272312, 0.06270845678671, 0.048590695510175]]

Example 4. Output Structure of a Multi-State Process

During model self-checks, you want to confirm that the variance paths are always nonnegative and that the state dimensions are arranged by column. For a multi-state process, each path slice is a numPaths × stateDim matrix; column 0 is price, and column 1 is variance.

def hestonEvolve(t, state, dt, noise, p) {
    variance = max(state[, 1], 0.0)
    nextSpot = state[, 0] * exp((p["rate"] - 0.5 * variance) * dt + sqrt(variance * dt) * noise[, 0])
    nextVariance = max(variance + p["kappa"] * (p["theta"] - variance) * dt + p["xi"] * sqrt(variance * dt) * noise[, 1], 0.0)
    return matrix(nextSpot, nextVariance)
}
p = dict(STRING, ANY)
p["rate"] = 0.03
p["kappa"] = 2.0
p["theta"] = 0.04
p["xi"] = 0.5
corrMatrix = matrix(1.0 0.0, 0.0 1.0)
corrMatrix[0, 1] = -0.7
corrMatrix[1, 0] = -0.7
options = dict(STRING, ANY)
options["seed"] = 42
options["batchSize"] = 5000
options["corr"] = corrMatrix
simulation = mcSimulate(hestonEvolve, 100.0 0.04, 0..12 / 12.0, 20000, 2, p, options)
spotPath = simulation["paths"][size(simulation["paths"]) - 1]
[size(simulation["paths"]), rows(spotPath), cols(spotPath), avg(spotPath[, 1])]
// output: [13, 20000, 2, 0.044645538557127]

Related Functions: monteCarlo, mcEstimate