Introduction to StochSimR

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Ayush Kundu

2026-05-28

Overview

StochSimR is a modular simulation engine for stochastic processes. It provides a unified interface for simulating, analysing, and visualising paths from a wide range of processes used in finance, physics, queueing theory, and applied probability.

Every simulator returns a stoch_path object that plugs directly into the plotting and analysis functions.

library(StochSimR)

Classic Processes

Brownian Motion

Standard Brownian motion (Wiener process) is the building block for most continuous-time models. The sim_brownian() function supports both exact simulation and Brownian bridge interpolation.

bm <- sim_brownian(T_max = 1, n_steps = 1000, mu = 0, sigma = 1,
                   n_paths = 100)
bm
#> StochSimR path object
#>   Process : Brownian Motion
#>   Method  : exact
#>   Paths   : 100
#>   Steps   : 1001
#>   Time    : [0, 1]
#>   Params  : mu = 0, sigma = 1, x0 = 0, T_max = 1
plot_paths(bm, max_paths = 50, show_mean = TRUE, show_bands = TRUE)

The red line is the sample mean across 100 paths; the shaded band shows the pointwise 2.5%–97.5% quantile envelope.

Poisson Process

The counting process for random arrivals. Supports both homogeneous (constant rate) and inhomogeneous (time-varying rate) specifications.

# Homogeneous
pp <- sim_poisson(T_max = 10, lambda = 3, n_paths = 20)
plot_paths(pp, show_bands = FALSE)

# Inhomogeneous with sinusoidal rate
pp_ih <- sim_poisson(
  T_max = 10,
  lambda = function(t) 3 + 2 * sin(2 * pi * t / 5),
  lambda_max = 5,
  n_paths = 20
)
plot_paths(pp_ih, show_bands = FALSE,
           title = "Inhomogeneous Poisson (sinusoidal rate)")

Markov Chain

Discrete-time Markov chains are simulated exactly by sampling from the transition matrix row at each step.

P <- matrix(c(0.7, 0.2, 0.1,
              0.3, 0.4, 0.3,
              0.2, 0.3, 0.5), nrow = 3, byrow = TRUE)

mc <- sim_markov(P, n_steps = 300, x0 = 1, n_paths = 10)
plot_paths(mc, show_bands = FALSE, show_mean = FALSE)

Verify convergence to the stationary distribution:

long_mc <- sim_markov(P, n_steps = 50000, x0 = 1, n_paths = 1)
freq <- table(factor(long_mc$values, levels = 1:3)) / length(long_mc$values)
freq
#> 
#>         1         2         3 
#> 0.4585308 0.2789344 0.2625347

Advanced Processes

Geometric Brownian Motion (GBM)

The standard model for stock prices. Two methods are available: "exact" uses the known log-normal solution; "euler" uses Euler–Maruyama discretisation.

stock <- sim_gbm(T_max = 1, n_steps = 252, mu = 0.08, sigma = 0.25,
                 x0 = 100, n_paths = 200)
plot_paths(stock, max_paths = 80)

plot_distribution(stock)

Ornstein–Uhlenbeck Process

A mean-reverting diffusion. With method = "exact", the known Gaussian transition density is used – no discretisation error.

ou <- sim_ou(T_max = 10, n_steps = 1000, theta = 2, mu = 5,
             sigma = 0.5, x0 = 0, n_paths = 50)
plot_paths(ou, show_mean = TRUE, show_bands = TRUE)

Note how all paths converge toward the long-run mean of 5.

Jump-Diffusion (Merton Model)

Combines GBM with compound Poisson jumps in log-price.

jd <- sim_jump_diffusion(
  T_max = 1, n_steps = 1000,
  mu = 0.05, sigma = 0.2,
  lambda = 3, mu_j = -0.02, sigma_j = 0.1,
  x0 = 100, n_paths = 50
)
plot_paths(jd, max_paths = 30)

Hawkes (Self-Exciting) Process

A point process where each event temporarily increases the rate of future events. The branching ratio alpha/beta controls the degree of clustering.

hk <- sim_hawkes(T_max = 50, mu = 0.5, alpha = 0.8, beta = 1.2,
                 n_paths = 15)
plot_paths(hk, show_bands = FALSE, show_mean = TRUE)

Levy Processes

Four families are available:

# Gamma subordinator (non-decreasing)
gam <- sim_levy(T_max = 5, n_steps = 500, type = "gamma",
                shape = 2, rate = 1, n_paths = 20)
plot_paths(gam, show_bands = FALSE)

# Variance-Gamma
vg <- sim_levy(T_max = 1, n_steps = 500, type = "variance_gamma",
               theta = -0.1, sigma = 0.3, nu = 0.5, n_paths = 30)
plot_paths(vg, max_paths = 20)

Summary Statistics

The path_summary() function computes terminal value distribution, extremes, increment moments, skewness, and excess kurtosis:

path_summary(stock)
#>            statistic        value
#> 1            n_paths 200.00000000
#> 2            n_steps 252.00000000
#> 3              T_max   1.00000000
#> 4      terminal_mean 109.07224277
#> 5        terminal_sd  28.53336114
#> 6    terminal_median 106.09773145
#> 7       terminal_q05  71.43366397
#> 8       terminal_q95 168.93447104
#> 9           max_mean 125.97579084
#> 10          min_mean  85.23655245
#> 11    increment_mean   0.03600096
#> 12      increment_sd   1.66925307
#> 13 terminal_skewness   0.75091748
#> 14 terminal_kurtosis   0.53570260

The plot_acf_paths() function shows the autocorrelation of path increments, useful for verifying independence (Brownian) or detecting dependence (OU):

plot_acf_paths(ou, path_idx = 1)

Variance Reduction

StochSimR includes four variance reduction techniques. Each returns the estimate, standard error, and the variance reduction factor versus crude Monte Carlo.

Control Variates

Use the terminal stock price (whose expectation is known analytically) as a control for pricing a European call:

h_call <- function(vals) pmax(vals[nrow(vals), ] - 100, 0)
g_term <- function(vals) vals[nrow(vals), ]
E_g <- 100 * exp(0.05)

cv <- vr_control_variate(sim_gbm, h = h_call, g = g_term, E_g = E_g,
  n_paths = 5000, T_max = 1, n_steps = 252,
  mu = 0.05, sigma = 0.2, x0 = 100)

cat("Estimate:", round(cv$estimate, 4), "\n")
#> Estimate: 10.8007
cat("SE:      ", round(cv$se, 4), "\n")
#> SE:       0.0822
cat("Reduction:", round(cv$reduction_factor, 1), "x\n")
#> Reduction: 6.7 x
cat("Correlation(h, g):", round(cv$correlation, 3), "\n")
#> Correlation(h, g): 0.922

Importance Sampling

Tilt the drift upward to estimate a tail probability more efficiently:

h_tail <- function(vals) as.numeric(vals[nrow(vals), ] > 150)
is_result <- vr_importance(sim_gbm, h = h_tail, drift_shift = 0.4,
  n_paths = 10000, T_max = 1, n_steps = 252,
  mu = 0.05, sigma = 0.2, x0 = 100)

cat("P(S_T > 150):", round(is_result$estimate, 6), "\n")
#> P(S_T > 150): 0.030385
cat("ESS:         ", round(is_result$ess), "/", is_result$n_samples, "\n")
#> ESS:          149 / 10000
cat("Reduction:   ", round(is_result$reduction_factor, 1), "x\n")
#> Reduction:    13.9 x

Monte Carlo Estimation

The mc_estimate() function provides a general-purpose estimator with batch-means standard errors and a convergence trace:

h_pos <- function(vals) pmax(vals[nrow(vals), ], 0)
mc <- mc_estimate(sim_brownian, h = h_pos, n_paths = 10000,
                  batch_size = 500, T_max = 1, n_steps = 200)

cat("E[max(B_1, 0)] =", round(mc$estimate, 4),
    "+/-", round(mc$se, 4), "\n")
#> E[max(B_1, 0)] = 0.3985 +/- 0.0058
cat("Exact answer:    ", round(1 / sqrt(2 * pi), 4), "\n")
#> Exact answer:     0.3989

Method Comparison

Compare simulation methods on the same process to measure bias and speed:

comp <- compare_methods(sim_ou, methods = c("exact", "euler"),
  n_paths = 500, n_reps = 20,
  T_max = 5, n_steps = 500, theta = 2, mu = 1, sigma = 0.5, x0 = 0)
print(comp)
#>   method mean_terminal sd_terminal bias rmse time_seconds
#> 1  exact     1.0003178  0.01585589   NA   NA      0.26365
#> 2  euler     0.9982826  0.01037739   NA   NA      0.22705

The exact method uses the known Gaussian transition density and has zero discretisation bias; the Euler method is faster per step but introduces \(O(\Delta t)\) weak error.

Session Info

sessionInfo()
#> R version 4.5.3 (2026-03-11)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Tahoe 26.5
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: Asia/Kolkata
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] StochSimR_1.1.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6        future.apply_1.20.2 jsonlite_2.0.0     
#>  [4] dplyr_1.2.0         compiler_4.5.3      tidyselect_1.2.1   
#>  [7] parallel_4.5.3      jquerylib_0.1.4     globals_0.19.1     
#> [10] scales_1.4.0        yaml_2.3.12         fastmap_1.2.0      
#> [13] ggplot2_4.0.3       R6_2.6.1            labeling_0.4.3     
#> [16] generics_0.1.4      knitr_1.51          future_1.70.0      
#> [19] tibble_3.3.1        bslib_0.10.0        pillar_1.11.1      
#> [22] RColorBrewer_1.1-3  rlang_1.2.0         cachem_1.1.0       
#> [25] xfun_0.57           sass_0.4.10         S7_0.2.1           
#> [28] otel_0.2.0          cli_3.6.6           withr_3.0.2        
#> [31] magrittr_2.0.4      digest_0.6.39       grid_4.5.3         
#> [34] rstudioapi_0.18.0   lifecycle_1.0.5     vctrs_0.7.1        
#> [37] evaluate_1.0.5      glue_1.8.0          farver_2.1.2       
#> [40] listenv_0.10.1      codetools_0.2-20    parallelly_1.46.1  
#> [43] rmarkdown_2.31      tools_4.5.3         pkgconfig_2.0.3    
#> [46] htmltools_0.5.9

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.