## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4,
  dpi = 100
)
set.seed(42)

## -----------------------------------------------------------------------------
library(StochSimR)

## ----brownian-----------------------------------------------------------------
bm <- sim_brownian(T_max = 1, n_steps = 1000, mu = 0, sigma = 1,
                   n_paths = 100)
bm
plot_paths(bm, max_paths = 50, show_mean = TRUE, show_bands = TRUE)

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

## ----poisson-inhom------------------------------------------------------------
# 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-------------------------------------------------------------------
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)

## ----markov-stationary--------------------------------------------------------
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

## ----gbm----------------------------------------------------------------------
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)

## ----ou-----------------------------------------------------------------------
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)

## ----jump-diffusion-----------------------------------------------------------
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-------------------------------------------------------------------
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-gamma---------------------------------------------------------------
# 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)

## ----levy-vg------------------------------------------------------------------
# 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------------------------------------------------------------------
path_summary(stock)

## ----acf----------------------------------------------------------------------
plot_acf_paths(ou, path_idx = 1)

## ----cv-----------------------------------------------------------------------
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")
cat("SE:      ", round(cv$se, 4), "\n")
cat("Reduction:", round(cv$reduction_factor, 1), "x\n")
cat("Correlation(h, g):", round(cv$correlation, 3), "\n")

## ----is-----------------------------------------------------------------------
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")
cat("ESS:         ", round(is_result$ess), "/", is_result$n_samples, "\n")
cat("Reduction:   ", round(is_result$reduction_factor, 1), "x\n")

## ----mc-----------------------------------------------------------------------
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")
cat("Exact answer:    ", round(1 / sqrt(2 * pi), 4), "\n")

## ----compare------------------------------------------------------------------
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)

## ----session------------------------------------------------------------------
sessionInfo()