mcEstimate

mcEstimate(samples, [confidenceLevel=0.95])

Computes the mean, unbiased sample variance, standard error of the mean, and normal asymptotic confidence interval for existing Monte Carlo samples.

samples is a numeric vector or numeric matrix, converted to DOUBLE type before computation. When the input is a numeric vector, it represents single-output samples; when it is a numeric matrix, it represents multi-output samples, where each row is a sample and each column is a target to be estimated. The number of samples must be at least 2, and the matrix must contain at least 1 column. All elements must be non-NULL finite numeric values.

confidenceLevel (optional) is a DOUBLE scalar that specifies the confidence level. The default is 0.95, and the valid range is (0,1).

Returns a dictionary containing the sample mean, unbiased sample variance, standard error of the mean, normal asymptotic confidence interval, and other results:

• When samples is a numeric vector, the first five fields are DOUBLE scalars; when samples is a numeric matrix, they are DOUBLE vectors with length equal to the number of input columns.
• The confidence interval is the estimate plus or minus the product of the corresponding standard normal quantile and the standard error, where the quantile is (1 + confidenceLevel) / 2.
• For multi-column samples, marginal confidence intervals are computed column by column; this does not represent a joint confidence interval that covers all outputs simultaneously.

Example 1. Compute Statistics for a Single-Output Sample

The first call computes statistics for 5 samples. The second call estimates the two-column sample column by column at a 99% confidence level. The results are, in order, the mean, unbiased sample variance, standard error, and sample count; for multi-column estimates, they also include the means and the lower and upper confidence interval bounds.

stats = mcEstimate(1.0 2.0 3.0 4.0 5.0)
samples = matrix(1.0 2.0 3.0 4.0, 10.0 20.0 30.0 40.0)
multi = mcEstimate(samples, 0.99)
[stats["estimate"], stats["variance"], stats["stdError"], stats["sampleCount"], multi["estimate"], multi["confidenceLower"], multi["confidenceUpper"]]
// output: [3, 2.5, 0.707106781186548, 5, [2.5, 25], [0.837309334113677, 8.37309334113677], [4.162690665886324, 41.626906658863234]]

Example 2. Compute Statistics Column by Column for Multiple Outputs

The input is a 4-by-2 matrix with a specified 99% confidence level. Returns the mean of the two columns.

samples = matrix(1.0 2.0 3.0 4.0, 10.0 20.0 30.0 40.0)
result = mcEstimate(samples=samples, confidenceLevel=0.99)
result["estimate"]
// output: [2.5, 25]

Example 3. Pair Antithetic Samples Before Computing Statistics

When antithetic=true is specified, two adjacent paths are antithetic to each other and form one statistical unit. If you pass all samples directly as independent samples, the estimate remains unchanged, but the standard error is overestimated and the confidence interval is too wide. The correct approach is to first average each pair, then pass the averages.

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["antithetic"] = true
options["batchSize"] = 5000
simulation = mcSimulate(gbmEvolve, 100.0, 0..12 / 12.0, 20000, 1, p, options)
terminal = flatten(simulation["paths"][size(simulation["paths"]) - 1])
samples = exp(-0.03) * max(terminal - 100.0, 0.0)
naive = mcEstimate(samples)
idx = 0..(size(samples) / 2 - 1)
pairedSamples = (samples[2 * idx] + samples[2 * idx + 1]) / 2.0
paired = mcEstimate(pairedSamples)
[naive["estimate"], paired["estimate"], naive["stdError"], paired["stdError"]]
// output: [9.350371878604676, 9.350371878604683, 0.099102747780388, 0.073832573115757]

Example 4. Interval Width for Different Confidence Levels

Internal risk control and regulatory reporting often use different conventions. The point estimate and standard error do not change with the confidence level; only the interval endpoints move.

samples = 1.2 0.8 1.5 1.1 0.9 1.3 1.0 1.4
wide = mcEstimate(samples, 0.99)
narrow = mcEstimate(samples, 0.90)
[wide["confidenceLower"], wide["confidenceUpper"], narrow["confidenceLower"], narrow["confidenceUpper"]]
// output: [0.926926638731427, 1.373073361268572, 1.007551497355305, 1.292448502644694]

Related Functions: monteCarlo, mcSimulate