## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)

## ----setup, echo=FALSE, message=FALSE-----------------------------------------
library(RandomWalker)
library(dplyr)
library(ggplot2)

## ----summary_basic------------------------------------------------------------
# Generate walks
walks <- random_normal_walk(.num_walks = 30, .n = 100)

# Overall summary
walks |> summarize_walks(.value = y)

## ----summary_by_walk----------------------------------------------------------
# Summary by walk
walks |>
  summarize_walks(.value = y, .group_var = walk_number) |>
  head()

## ----summary_different_values-------------------------------------------------
# Summarize cumulative sum
walks |> summarize_walks(.value = cum_sum_y)

## ----summary_gbm--------------------------------------------------------------
# Summarize cumulative product
geometric_brownian_motion(.num_walks = 30, .initial_value = 100) |>
  summarize_walks(.value = cum_prod_y)

## ----summary_by_group---------------------------------------------------------
# Summarize by group
walks |>
  summarize_walks(.value = cum_sum_y, .group_var = walk_number) |>
  head()

## ----stock_simulation---------------------------------------------------------
# Simulate stock prices
stock_sim <- geometric_brownian_motion(
  .num_walks = 1000,
  .n = 252,  # Trading days
  .mu = 0.08,
  .sigma = 0.25,
  .initial_value = 100
)

# Get final price statistics
final_prices <- stock_sim |>
  summarize_walks(.value = cum_prod_y, .group_var = walk_number) |>
  pull(max_val)

# Analyze outcomes
tibble(final_price = final_prices) |>
  summarize(
    median_price = median(final_price),
    mean_price = mean(final_price),
    prob_profit = mean(final_price > 100),
    prob_loss_20 = mean(final_price < 80),
    sd_returns = sd((final_price - 100) / 100)
  )

## ----compare_distributions----------------------------------------------------
# Normal vs Cauchy walks
normal_stats <- random_normal_walk(.num_walks = 100, .n = 100) |>
  summarize_walks(.value = y) |>
  mutate(distribution = "Normal")

cauchy_stats <- random_cauchy_walk(.num_walks = 100, .n = 100) |>
  summarize_walks(.value = y) |>
  mutate(distribution = "Cauchy")

# Compare
bind_rows(normal_stats, cauchy_stats) |>
  select(distribution, mean_val, sd, skewness, kurtosis)

## ----cumulative_functions-----------------------------------------------------
# Generate walk
walks <- random_normal_walk(.num_walks = 10, .n = 100, .initial_value = 100)

# Cumulative functions are already in the data
walks |>
  select(walk_number, step_number, y, starts_with("cum_")) |>
  head(10)

## ----analyze_cumsum-----------------------------------------------------------
# Analyze cumulative sum
walks |>
  summarize_walks(.value = cum_sum_y, .group_var = walk_number) |>
  head()

## ----track_extremes-----------------------------------------------------------
# Track maximum ever reached
walks |>
  group_by(walk_number) |>
  summarize(
    max_ever = max(cum_max_y),
    min_ever = min(cum_min_y),
    final_value = last(cum_sum_y)
  ) |>
  head()

## ----custom_cumulative, fig.alt="Line plot showing maximum drawdown over time for 10 random walks. Each walk is shown in a different color with semi-transparent lines. The y-axis shows the maximum drawdown as a negative percentage, and the x-axis shows the step number from 0 to 100."----
# Add custom cumulative functions
walks_extended <- walks |>
  group_by(walk_number) |>
  mutate(
    # Cumulative variance
    cum_var = cumsum((y - cumsum(y) / row_number())^2) / row_number(),
    # Cumulative absolute sum
    cum_abs_sum = cumsum(abs(y)),
    # Running maximum drawdown
    running_peak = cummax(cum_sum_y),
    drawdown = (cum_sum_y - running_peak) / running_peak,
    max_drawdown = cummin(drawdown)
  ) |>
  ungroup()

# Visualize drawdown
walks_extended |>
  ggplot(aes(x = step_number, y = max_drawdown, color = walk_number)) +
  geom_line(alpha = 0.7) +
  theme_minimal() +
  labs(title = "Maximum Drawdown Over Time")

## ----confidence_interval------------------------------------------------------
# Generate data
x <- rnorm(1000, mean = 10, sd = 2)

# Calculate 95% CI (default)
confidence_interval(x)

## ----confidence_interval_99---------------------------------------------------
# Calculate 99% CI
confidence_interval(x, .interval = 0.01)

## ----confidence_interval_90---------------------------------------------------
# Calculate 90% CI
confidence_interval(x, .interval = 0.10)

## ----ci_by_step, fig.alt="Time series plot showing the mean random walk trajectory (dark blue line) surrounded by a light blue shaded region representing the 95% confidence interval. The x-axis shows step numbers from 0 to 100, and the y-axis shows the walk values."----
# Generate walks
walks <- random_normal_walk(.num_walks = 100, .n = 100)

# Calculate CI at each step
ci_by_step <- walks |>
  group_by(step_number) |>
  summarize(
    mean_val = mean(y),
    lower = quantile(y, 0.025),
    upper = quantile(y, 0.975)
  )

# Visualize
ggplot(ci_by_step, aes(x = step_number)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.3, fill = "steelblue") +
  geom_line(aes(y = mean_val), color = "darkblue", linewidth = 1) +
  theme_minimal() +
  labs(
    title = "Mean Random Walk with 95% Confidence Interval",
    x = "Step",
    y = "Value"
  )

## ----ci_final_values----------------------------------------------------------
# Get final values from many walks
walks <- random_normal_walk(.num_walks = 1000, .n = 100, .initial_value = 100)

final_values <- walks |>
  group_by(walk_number) |>
  slice_max(step_number, n = 1) |>
  pull(cum_sum_y)

# Calculate confidence interval
confidence_interval(final_values)

## ----running_quantile---------------------------------------------------------
# Generate walks
walks <- random_normal_walk(.num_walks = 100, .n = 100)

# Calculate running median (50th percentile)
walks_with_median <- walks |>
  group_by(step_number) |>
  mutate(median_at_step = running_quantile(y, .probs = 0.5, .window = 5)) |>
  ungroup()

# Show results
walks_with_median |>
  select(walk_number, step_number, y, median_at_step) |>
  head(10)

## ----running_quartiles--------------------------------------------------------
# Calculate running quartiles
walks_with_quartiles <- walks |>
  group_by(step_number) |>
  mutate(
    q25 = running_quantile(y, .probs = 0.25, .window = 5),
    q50 = running_quantile(y, .probs = 0.50, .window = 5),
    q75 = running_quantile(y, .probs = 0.75, .window = 5)
  ) |>
  ungroup()

# Show results
walks_with_quartiles |>
  select(walk_number, step_number, y, q25, q50, q75) |>
  head(10)

## ----quantile_evolution, fig.alt="Plot showing the evolution of quantiles over 100 steps for 200 random walks. A dark blue line shows the median (50th percentile), surrounded by darker blue shading showing the interquartile range (25th to 75th percentile), and lighter blue shading showing the 90% confidence interval (5th to 95th percentile). The x-axis shows step numbers and the y-axis shows walk values."----
# Generate many walks
walks <- random_normal_walk(.num_walks = 200, .n = 100)

# Calculate quantiles at each step
quantile_evolution <- walks |>
  group_by(step_number) |>
  summarize(
    q05 = quantile(y, 0.05),
    q25 = quantile(y, 0.25),
    q50 = quantile(y, 0.50),
    q75 = quantile(y, 0.75),
    q95 = quantile(y, 0.95)
  )

# Plot
ggplot(quantile_evolution, aes(x = step_number)) +
  geom_ribbon(aes(ymin = q05, ymax = q95), alpha = 0.2, fill = "blue") +
  geom_ribbon(aes(ymin = q25, ymax = q75), alpha = 0.3, fill = "blue") +
  geom_line(aes(y = q50), color = "darkblue", linewidth = 1) +
  theme_minimal() +
  labs(
    title = "Random Walk Quantile Evolution",
    subtitle = "Median (dark blue), IQR (darker shading), and 90% CI (lighter shading)",
    x = "Step",
    y = "Value"
  )

## ----euclidean_distance, fig.alt="Multi-line plot showing Euclidean distance from origin over time for 10 different 2D random walks. Each walk is shown in a different color with semi-transparent lines. The x-axis shows step numbers from 0 to 100, and the y-axis shows the distance from the origin."----
# 2D walk
walks_2d <- random_normal_walk(.num_walks = 10, .n = 100, .dimensions = 2)

# Calculate Euclidean distance
walks_with_distance <- walks_2d |>
  euclidean_distance(.x = x, .y = y)

# Visualize distance over time
walks_with_distance |>
  ggplot(aes(x = step_number, y = distance, color = walk_number)) +
  geom_line(alpha = 0.7) +
  theme_minimal() +
  labs(
    title = "Distance from Origin in 2D Random Walk",
    x = "Step",
    y = "Euclidean Distance"
  )

## ----distance_stats, fig.alt="Line plot comparing observed mean distance (solid blue line) versus theoretical prediction (dashed red line) for 3D random walks. The theoretical prediction follows sqrt(n) as expected for standard Brownian motion. The x-axis shows step numbers and the y-axis shows distance from origin."----
# 3D walk
walks_3d <- random_normal_walk(.num_walks = 100, .n = 1000, .dimensions = 3)

# Calculate distance
walks_with_dist <- walks_3d |> euclidean_distance(.x = x, .y = z)

# Analyze distance evolution
distance_stats <- walks_with_dist |>
  group_by(step_number) |>
  summarize(
    mean_dist = mean(distance),
    sd_dist = sd(distance),
    max_dist = max(distance)
  )

# Plot average distance vs sqrt(n) theoretical prediction
distance_stats |>
  ggplot(aes(x = step_number)) +
  geom_line(aes(y = mean_dist, color = "Observed"), linewidth = 1) +
  geom_line(aes(y = sqrt(3 * step_number), color = "Theory"), linewidth = 1, linetype = "dashed") +
  scale_color_manual(values = c("Observed" = "blue", "Theory" = "red")) +
  theme_minimal() +
  labs(
    title = "Mean Distance vs Theoretical Prediction",
    subtitle = "Distance ~ sqrt(3n) for 3D Brownian motion",
    x = "Step",
    y = "Distance",
    color = ""
  )

## ----first_passage, fig.alt="Histogram showing the distribution of first passage times to level 10 for 100 discrete random walks. The x-axis shows the first passage time in steps, and the y-axis shows the count of walks. The distribution is right-skewed with most walks reaching level 10 within the first few hundred steps."----
# Generate walks
walks <- discrete_walk(.num_walks = 100, .n = 1000, .initial_value = 0)

# Find first passage time to level 10
first_passage <- walks |>
  group_by(walk_number) |>
  filter(cum_sum_y >= 10) |>
  slice_min(step_number, n = 1) |>
  select(walk_number, first_passage_time = step_number)

# Analyze distribution of first passage times
first_passage |>
  ggplot(aes(x = first_passage_time)) +
  geom_histogram(bins = 50, fill = "steelblue", alpha = 0.7) +
  theme_minimal() +
  labs(
    title = "Distribution of First Passage Times to Level 10",
    x = "First Passage Time (Steps)",
    y = "Count"
  )

## ----subset_walks-------------------------------------------------------------
# Generate walks
walks <- random_normal_walk(.num_walks = 100, .n = 100, .initial_value = 100)

# Get walk with maximum final value
max_walk <- walks |> subset_walks(.value = "cum_sum_y", .type = "max")

# Get walk with minimum final value
min_walk <- walks |> subset_walks(.value = "cum_sum_y", .type = "min")

# Show the extreme walks
max_walk |> 
  summarize_walks(.value = cum_sum_y, .group_var = walk_number)

min_walk |> 
  summarize_walks(.value = cum_sum_y, .group_var = walk_number)

## ----find_specific_walks, fig.alt="Multi-line plot showing all random walks that crossed the threshold of 102 at some point during their evolution. Multiple semi-transparent lines in different colors show the cumulative sum trajectories. The x-axis shows step numbers from 0 to 100, and the y-axis shows the cumulative sum values."----
# Find walks that cross a threshold
walks <- random_normal_walk(.num_walks = 100, .n = 100, .initial_value = 100)

# Identify walks that reached 102
crossed_102 <- walks |>
  group_by(walk_number) |>
  filter(any(cum_sum_y >= 102)) |>
  pull(walk_number) |>
  unique()

# Extract and visualize those walks
walks |>
  filter(walk_number %in% crossed_102) |>
  visualize_walks(.pluck = "cum_sum_y", .alpha = 0.3)

## ----autocorrelation, fig.alt="Autocorrelation function (ACF) plot showing the correlation of random walk steps with their lagged values. Vertical bars represent correlation values at different lags, with dashed blue lines indicating significance thresholds."----
# Generate walk with drift
walks <- random_normal_drift_walk(.num_walks = 1, .n = 500, .mu = 0.1)

# Calculate autocorrelation
acf_result <- walks |> pull(y) |> acf(plot = FALSE)

# Plot
plot(acf_result, main = "Autocorrelation of Random Walk Steps")

## ----distribution_testing, fig.alt="Q-Q (quantile-quantile) plot comparing the distribution of random walk steps to a theoretical normal distribution. Points falling along the diagonal reference line indicate good agreement with normality."----
# Generate walks
walks <- random_normal_walk(.num_walks = 100, .n = 100)

# Test if steps are normally distributed
steps <- walks |> pull(y)

# Shapiro-Wilk test for normality
shapiro.test(sample(steps, 5000))  # Sample for computational efficiency

# Q-Q plot
qqnorm(steps)
qqline(steps)

## ----variance_ratio-----------------------------------------------------------
# Generate walk
walk <- random_normal_walk(.num_walks = 1, .n = 1000)

# Calculate variance ratio
values <- walk |> pull(cum_sum_y)

# Variance of k-differences
k <- 10
var_k <- var(diff(values, lag = k))
var_1 <- var(diff(values, lag = 1))

# Variance ratio (should be ≈ k for random walk)
vr <- var_k / (k * var_1)
print(paste("Variance Ratio:", round(vr, 3), "| Expected:", k))

## ----return_analysis----------------------------------------------------------
# Generate stock price simulation
prices <- geometric_brownian_motion(
  .num_walks = 1,
  .n = 252,
  .mu = 0.08,
  .sigma = 0.25,
  .initial_value = 100
)

# Calculate returns
returns <- prices |>
  mutate(
    log_return = log(cum_prod_y / lag(cum_prod_y)),
    simple_return = (cum_prod_y - lag(cum_prod_y)) / lag(cum_prod_y)
  ) |>
  filter(!is.na(log_return))

# Analyze returns
returns |>
  summarize(
    mean_return = mean(log_return) * 252,  # Annualized
    volatility = sd(log_return) * sqrt(252),  # Annualized
    sharpe_ratio = mean_return / volatility
  )

## ----compare_tests------------------------------------------------------------
# Generate two types of walks
normal_walks <- random_normal_walk(.num_walks = 50, .n = 100)
cauchy_walks <- random_cauchy_walk(.num_walks = 50, .n = 100)

# Get final values
normal_final <- normal_walks |>
  group_by(walk_number) |>
  slice_max(step_number) |>
  pull(cum_sum_y)

cauchy_final <- cauchy_walks |>
  group_by(walk_number) |>
  slice_max(step_number) |>
  pull(cum_sum_y)

# Wilcoxon rank-sum test (non-parametric)
wilcox.test(normal_final, cauchy_final)

## ----ks_test------------------------------------------------------------------
# Kolmogorov-Smirnov test
ks.test(normal_final, cauchy_final)

## ----test_drift---------------------------------------------------------------
# Generate walk with known drift
walks <- random_normal_drift_walk(.num_walks = 100, .n = 100, .mu = 0.1)

# Test if mean step is significantly different from 0
steps <- walks |> pull(y)
t.test(steps, mu = 0)