Simulate Stochastic Process

Simulate sample paths from stochastic processes -- including discrete Markov chains, continuous-time processes, stochastic differential equations, and MCMC samplers -- with convergence diagnostics, variance reduction techniques, and trajectory visualization.

When to Use

• You need to generate sample paths from a stochastic process for estimation, prediction, or visualization
• Analytical solutions are intractable and simulation is the only feasible approach
• You are running Monte Carlo estimation and need convergence guarantees and uncertainty quantification
• You want to validate analytical results (stationary distributions, hitting times) against empirical simulation
• You need to sample from a complex posterior distribution using MCMC
• You are prototyping a stochastic model before committing to full analytical treatment

Inputs

Required

process_type

"dtmc"

"ctmc"

"random_walk"

"brownian_motion"

"sde"

"mcmc"

parameters

n_paths

n_steps

Optional

initial_state

dt

seed

burn_in

n_steps / 10

thinning

variance_reduction

"none"

"none"

"antithetic"

"stratified"

"control_variate"

target_function

Procedure

Step 1: Define Process Model and Parameters

1.1. Identify the process type and gather all required parameters:

• DTMC: Transition matrix

  P

  and state space. Validate

  P

  is row-stochastic.
• CTMC: Rate matrix

  Q

  . Validate rows sum to 0 and off-diagonal entries are non-negative.
• Random walk: Step distribution (e.g.,

  {-1, +1}

  with equal probability), boundaries if any.
• Brownian motion: Drift

  mu

  , volatility

  sigma

  , dimension

  d

  .
• SDE (Ito): Drift function

  a(x,t)

  , diffusion function

  b(x,t)

  .
• MCMC: Target log-density, proposal mechanism (random walk Metropolis, Hamiltonian, Gibbs components).

1.2. Validate parameter consistency:

• Matrix dimensions match state space size.
• SDE coefficients satisfy growth and Lipschitz conditions (at least informally) for the chosen solver.
• MCMC proposal is well-defined for the support of the target distribution.

1.3. Set the random seed for reproducibility.

Expected: A fully specified stochastic model with validated parameters and a reproducible random state.

On failure: If parameters are inconsistent (e.g., non-stochastic matrix), correct them before proceeding. If SDE coefficients are pathological, consider a different discretization scheme.

Step 2: Select Simulation Method

2.1. Choose the appropriate algorithm based on process type:

dt

2.2. For SDE methods, choose

dt

dt = 0.01

2.3. For MCMC, tune the proposal scale to achieve an acceptance rate of approximately:

• 23.4% for high-dimensional random walk Metropolis
• 57.4% for one-dimensional targets
• 65-90% for HMC (depends on trajectory length)

2.4. If variance reduction is requested, configure it:

• Antithetic variates: For each path with random increments

  Z

  , also simulate with

  -Z

  .
• Stratified sampling: Partition the probability space and sample within each stratum.
• Control variates: Identify a correlated quantity with known expectation to reduce variance.

Expected: A selected simulation algorithm matched to the process type with appropriate tuning parameters.

On failure: If the chosen method is unstable (e.g., Euler-Maruyama diverging), switch to an implicit method or reduce

dt

Step 3: Implement and Run Simulation

3.1. Allocate storage for

n_paths

n_steps

3.2. For each path

i = 1, ..., n_paths

DTMC / Random Walk:

• Set

  x[0] = initial_state

• For

  t = 1, ..., n_steps

  : sample

  x[t]

  from the transition distribution given

  x[t-1]

CTMC (Gillespie):

• Set

  x[0] = initial_state

  ,

  time = 0

• While

  time < T_max

  :
  • Compute total rate

    lambda = -Q[x, x]

  • Sample holding time

    tau ~ Exp(lambda)

  • Sample next state from transition probabilities

    Q[x, j] / lambda

    for

    j != x

  • Update

    time += tau

    , record transition

SDE (Euler-Maruyama):

• Set

  x[0] = initial_state

• For

  t = 1, ..., n_steps

  :

  • dW = sqrt(dt) * N(0, I)

    (Wiener increment)

  • x[t] = x[t-1] + a(x[t-1], t*dt) * dt + b(x[t-1], t*dt) * dW

MCMC (Metropolis-Hastings):

• Set

  x[0] = initial_state

• For

  t = 1, ..., n_steps

  :
  • Propose

    x' ~ q(x' | x[t-1])

  • Compute acceptance ratio

    alpha = min(1, p(x') * q(x[t-1]|x') / (p(x[t-1]) * q(x'|x[t-1])))

  • Accept with probability

    alpha

    :

    x[t] = x'

    if accepted, else

    x[t] = x[t-1]

  • Record acceptance decision

3.3. If

target_function

3.4. Apply thinning: keep every

thinning

3.5. Discard

burn_in

Expected:

n_paths

On failure: If simulation produces NaN or Inf values, reduce

dt

Step 4: Apply Convergence Diagnostics

4.1. Trace plots: Plot the value of each component over time for a subset of paths. Visual inspection for stationarity (no trends, stable variance).

4.2. Gelman-Rubin diagnostic (R-hat): For MCMC with multiple chains:

• Compute within-chain variance

  W

  and between-chain variance

  B

  .

• R_hat = sqrt((n-1)/n + B/(n*W))

• Convergence indicated by

  R_hat < 1.01

  (strict) or

  R_hat < 1.1

  (lenient).

4.3. Effective sample size (ESS):

• Estimate autocorrelation at increasing lags.

• ESS = n_samples / (1 + 2 * sum(autocorrelations))

• Rule of thumb:

  ESS > 400

  for reliable posterior summaries.

4.4. Geweke diagnostic: Compare the mean of the first 10% and last 50% of each chain. The z-score should be within [-2, 2] for convergence.

4.5. For non-MCMC processes: Verify that time-averaged statistics (mean, variance) stabilize as path length increases. Plot running averages.

4.6. Report a summary table:

Expected: All convergence diagnostics pass their thresholds. Trace plots show stable, well-mixing chains.

On failure: If R-hat > 1.1, run longer chains or improve the proposal. If ESS is very low, increase thinning or switch to a better sampler (e.g., HMC). If Geweke fails, extend burn-in.

Step 5: Compute Summary Statistics with Confidence Intervals

5.1. For each quantity of interest (state occupancy, function expectation, hitting times):

• Compute the point estimate as the sample mean across paths (after burn-in and thinning).
• Compute the standard error using the effective sample size:

  SE = SD / sqrt(ESS)

  .

5.2. Construct confidence intervals:

• Normal approximation:

  estimate +/- z_{alpha/2} * SE

• For skewed distributions, use percentile bootstrap or batch means.

5.3. If variance reduction was applied, compute the variance reduction factor:

• VRF = Var(naive estimator) / Var(reduced estimator)

• Report the effective speedup.

5.4. For Monte Carlo integration estimates:

• Report the estimate, standard error, 95% CI, ESS, and number of function evaluations.

5.5. For distribution estimates:

• Compute empirical quantiles (median, 2.5th, 97.5th percentiles).
• Kernel density estimates for continuous quantities.

5.6. Tabulate all summary statistics with their uncertainties.

Expected: Point estimates with associated standard errors and confidence intervals. Variance reduction (if applied) yields a VRF > 1.

On failure: If confidence intervals are too wide, increase

n_paths

n_steps

Step 6: Visualize Trajectories and Distributions

6.1. Trajectory plots: Plot a representative subset of sample paths (5-20 paths) over time. Use transparency for overlapping paths.

6.2. Ensemble statistics: Overlay the mean trajectory and pointwise 95% confidence bands across all paths.

6.3. Marginal distributions: At selected time points, plot histograms or density estimates of the state distribution across paths.

6.4. Stationary distribution comparison: If an analytical stationary distribution is available, overlay it on the empirical histogram from the final time slice.

6.5. Autocorrelation plots: For MCMC, plot the autocorrelation function (ACF) for each component up to a reasonable lag.

6.6. Diagnostic dashboard: Combine trace plots, ACF plots, running mean plots, and marginal densities into a single multi-panel figure for comprehensive assessment.

6.7. Save all figures in both vector (PDF/SVG) and raster (PNG) formats for documentation.

Expected: Publication-quality figures showing trajectory behavior, distributional convergence, and diagnostic summaries. Analytical solutions (where available) match empirical results.

On failure: If visualizations reveal non-stationarity or multimodality not expected from the model, revisit Steps 1-2 for parameter or method errors. If plots are cluttered, reduce the number of displayed paths or increase figure size.

Validation

• All simulated trajectories remain in the valid state space (no out-of-bounds values, no NaN/Inf)
• For DTMC/CTMC: empirical stationary distribution converges to the analytical one (within expected Monte Carlo error)
• For SDE: halving

  dt

  does not qualitatively change the results (convergence order check)
• For MCMC: R-hat < 1.01, ESS > 400, Geweke z-scores within [-2, 2]
• Confidence interval widths decrease proportionally to

  1/sqrt(n_paths)

  (central limit theorem)
• Variance reduction techniques yield VRF > 1 (estimates improve, not worsen)
• Reproducibility: re-running with the same seed produces identical results

Common Pitfalls

• Insufficient burn-in for MCMC: Starting from a poor initial state requires a long burn-in before samples represent the target distribution. Always inspect trace plots and use convergence diagnostics rather than guessing the burn-in length.
• Euler-Maruyama instability for stiff SDEs: If the drift term has large gradients, explicit Euler-Maruyama can diverge. Switch to implicit methods or use adaptive step sizing.
• Confusing strong and weak convergence for SDEs: Strong convergence measures pathwise error (important for individual trajectories); weak convergence measures distributional error (sufficient for expectations). Euler-Maruyama has weak order 1.0 but strong order 0.5.
• Pseudorandom number generator quality: For very long simulations, low-quality RNGs can produce correlated samples. Use a well-tested generator (Mersenne Twister, PCG, or Xoshiro) and verify independence.
• Ignoring autocorrelation in MCMC: Treating autocorrelated MCMC samples as independent underestimates uncertainty. Always use effective sample size, not raw sample count, for standard errors.
• Antithetic variates for non-monotone functions: Antithetic sampling reduces variance only when the estimand is a monotone function of the underlying uniforms. For non-monotone functions, it can increase variance.
• Memory for large simulations: Storing all time steps of many long paths can exhaust memory. Use online statistics (running mean, variance) when full trajectories are not needed for visualization.

Related Skills

• Model Markov Chain -- provides the transition matrices and analytical solutions that simulation validates
• Fit Hidden Markov Model -- simulation from fitted HMMs enables posterior predictive checking and synthetic data generation