library(ggplot2)
library(pracma) # For Hankel matrix construction
# Function to generate Fractional Brownian Motion via Cholesky decomposition
generate_fBM <- function(n, H) {
# Create covariance matrix
cov_matrix <- matrix(0, n, n)
for (i in 1:n) {
for (j in 1:n) {
cov_matrix[i, j] <- 0.5 * (abs(i)^(2*H) + abs(j)^(2*H) - abs(i-j)^(2*H))
}
}
# Cholesky decomposition
L <- chol(cov_matrix)
# Generate standard normal random variables
Z <- rnorm(n)
# Compute fBM
fBM <- L %*% Z
return(as.vector(fBM))
}
# Simulation parameters
n_samples <- 1000 # Number of samples
sample_size <- 100 # Length of each fBM path
hurst_values <- c(0.5, 0.6, 0.7, 0.8) # Hurst parameters to test
# Initialize results data frame
results <- data.frame(Hurst = numeric(), Variance = numeric())
# Perform simulations
set.seed(123) # For reproducibility
for (H in hurst_values) {
for (i in 1:n_samples) {
# Generate fBM sample
fbm_sample <- generate_fBM(sample_size, H)
# Estimate variance (standard approach for Brownian motion)
var_estimate <- var(diff(fbm_sample))
# Store results
results <- rbind(results, data.frame(Hurst = H, Variance = var_estimate))
}
}
# Theoretical variances (for standard fBM increments)
theoretical_vars <- c(1, 1, 1, 1) # Since we're looking at increments
# Create boxplots with light blue color
ggplot(results, aes(x = factor(Hurst), y = Variance, fill = factor(Hurst))) +
geom_boxplot(fill = "lightblue", color = "blue", alpha = 0.7) +
geom_hline(yintercept = 1, linetype = "dashed", color = "red", linewidth = 0.8) +
labs(title = "Variance Estimation of Fractional Brownian Motion Increments",
subtitle = paste("Based on", n_samples, "samples of size", sample_size),
x = "Hurst Parameter (H)",
y = "Estimated Variance",
caption = "Dashed red line shows theoretical variance of increments") +
scale_x_discrete(labels = paste0("H = ", hurst_values)) +
theme_minimal() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5),
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_blank())
library(ggplot2)
library(gridExtra) # For arranging multiple plots
library(grid) # For textGrob and gpar functions
# Function to generate Fractional Brownian Motion via Cholesky decomposition
generate_fBM <- function(n, H) {
# Create covariance matrix
cov_matrix <- matrix(0, n, n)
for (i in 1:n) {
for (j in 1:n) {
cov_matrix[i, j] <- 0.5 * (abs(i)^(2*H) + abs(j)^(2*H) - abs(i-j)^(2*H))
}
}
# Cholesky decomposition with error handling
L <- tryCatch(
chol(cov_matrix),
error = function(e) {
# Add small noise to ensure positive definiteness
cov_matrix <- cov_matrix + diag(n) * 1e-10
chol(cov_matrix)
}
)
# Generate standard normal random variables
Z <- rnorm(n)
# Compute fBM
fBM <- L %*% Z
return(as.vector(fBM))
}
# Simulation parameters
n_samples <- 1000 # Number of samples
sample_size <- 100 # Length of each fBM path
hurst_group1 <- c(0.5, 0.6, 0.7)
hurst_group2 <- c(0.8, 0.9)
# Initialize results data frames
results_group1 <- data.frame(Hurst = numeric(), Variance = numeric())
results_group2 <- data.frame(Hurst = numeric(), Variance = numeric())
# Perform simulations
set.seed(123) # For reproducibility
# Simulation for group 1 (H = 0.5, 0.6, 0.7)
for (H in hurst_group1) {
for (i in 1:n_samples) {
fbm_sample <- generate_fBM(sample_size, H)
var_estimate <- var(diff(fbm_sample)) # Variance of increments
results_group1 <- rbind(results_group1,
data.frame(Hurst = H, Variance = var_estimate))
}
}
# Simulation for group 2 (H = 0.8, 0.9)
for (H in hurst_group2) {
for (i in 1:n_samples) {
fbm_sample <- generate_fBM(sample_size, H)
var_estimate <- var(diff(fbm_sample)) # Variance of increments
results_group2 <- rbind(results_group2,
data.frame(Hurst = H, Variance = var_estimate))
}
}
# Create color palettes
blue_palette <- c("#E6F0FF", "#B3D1FF", "#80B3FF") # Light to medium blues
orange_palette <- c("#FFE6CC", "#FFB366", "#FF8000") # Light to medium oranges
# Create first plot (H = 0.5, 0.6, 0.7)
plot1 <- ggplot(results_group1, aes(x = factor(Hurst), y = Variance, fill = factor(Hurst))) +
geom_boxplot(alpha = 0.8, outlier.color = "blue", outlier.alpha = 0.3) +
geom_hline(yintercept = 1, linetype = "dashed", color = "red", linewidth = 0.8) +
scale_fill_manual(values = blue_palette) +
labs(title = "Variance Estimation for H = 0.5, 0.6, 0.7",
x = "Hurst Parameter (H)",
y = "Estimated Variance") +
theme_minimal() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.background = element_rect(fill = "white", color = NA),
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_blank())
# Create second plot (H = 0.8, 0.9)
plot2 <- ggplot(results_group2, aes(x = factor(Hurst), y = Variance, fill = factor(Hurst))) +
geom_boxplot(alpha = 0.8, outlier.color = "darkorange", outlier.alpha = 0.3) +
geom_hline(yintercept = 1, linetype = "dashed", color = "red", linewidth = 0.8) +
scale_fill_manual(values = orange_palette) +
labs(title = "Variance Estimation for H = 0.8, 0.9",
x = "Hurst Parameter (H)",
y = "Estimated Variance") +
theme_minimal() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.background = element_rect(fill = "white", color = NA),
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_blank())
# Arrange plots side by side with main title
grid.arrange(plot1, plot2, ncol = 2,
top = textGrob("Fractional Brownian Motion Variance Estimation\n(1000 samples, n=100 each)",
gp = gpar(fontsize = 14, fontface = "bold")))
library(ggplot2)
library(gridExtra) # For arranging multiple plots
library(grid) # For textGrob and gpar functions
# Function to generate Fractional Brownian Motion via Cholesky decomposition
generate_fBM <- function(n, H) {
# Create covariance matrix
cov_matrix <- matrix(0, n, n)
for (i in 1:n) {
for (j in 1:n) {
cov_matrix[i, j] <- 0.5 * (abs(i)^(2*H) + abs(j)^(2*H) - abs(i-j)^(2*H))
}
}
# Cholesky decomposition with error handling
L <- tryCatch(
chol(cov_matrix),
error = function(e) {
# Add small noise to ensure positive definiteness
cov_matrix <- cov_matrix + diag(n) * 1e-10
chol(cov_matrix)
}
)
# Generate standard normal random variables
Z <- rnorm(n)
# Compute fBM
fBM <- L %*% Z
return(as.vector(fBM))
}
# Simulation parameters
n_samples <- 200 # Number of samples
sample_size <- 100 # Length of each fBM path
hurst_group1 <- c(0.5, 0.6, 0.7)
hurst_group2 <- c(0.8, 0.9)
# Initialize results data frames
results_group1 <- data.frame(Hurst = numeric(), Variance = numeric())
results_group2 <- data.frame(Hurst = numeric(), Variance = numeric())
# Perform simulations
set.seed(123) # For reproducibility
# Simulation for group 1 (H = 0.5, 0.6, 0.7)
for (H in hurst_group1) {
for (i in 1:n_samples) {
fbm_sample <- generate_fBM(sample_size, H)
var_estimate <- var(diff(fbm_sample)) # Variance of increments
results_group1 <- rbind(results_group1,
data.frame(Hurst = H, Variance = var_estimate))
}
}
# Simulation for group 2 (H = 0.8, 0.9)
for (H in hurst_group2) {
for (i in 1:n_samples) {
fbm_sample <- generate_fBM(sample_size, H)
var_estimate <- var(diff(fbm_sample)) # Variance of increments
results_group2 <- rbind(results_group2,
data.frame(Hurst = H, Variance = var_estimate))
}
}
# Create color palettes
blue_palette <- c("#E6F0FF", "#B3D1FF", "#80B3FF") # Light to medium blues
orange_palette <- c("#FFE6CC", "#FFB366", "#FF8000") # Light to medium oranges
# Create first plot (H = 0.5, 0.6, 0.7)
plot1 <- ggplot(results_group1, aes(x = factor(Hurst), y = Variance, fill = factor(Hurst))) +
geom_boxplot(alpha = 0.8, outlier.color = "blue", outlier.alpha = 0.3) +
geom_hline(yintercept = 1, linetype = "dashed", color = "red", linewidth = 0.8) +
scale_fill_manual(values = blue_palette) +
labs(title = "Variance Estimation for H = 0.5, 0.6, 0.7",
x = "Hurst Parameter (H)",
y = "Estimated Variance") +
theme_minimal() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.background = element_rect(fill = "white", color = NA),
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_blank())
# Create second plot (H = 0.8, 0.9)
plot2 <- ggplot(results_group2, aes(x = factor(Hurst), y = Variance, fill = factor(Hurst))) +
geom_boxplot(alpha = 0.8, outlier.color = "darkorange", outlier.alpha = 0.3) +
geom_hline(yintercept = 1, linetype = "dashed", color = "red", linewidth = 0.8) +
scale_fill_manual(values = orange_palette) +
labs(title = "Variance Estimation for H = 0.8, 0.9",
x = "Hurst Parameter (H)",
y = "Estimated Variance") +
theme_minimal() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.background = element_rect(fill = "white", color = NA),
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_blank())
# Arrange plots side by side with main title
grid.arrange(plot1, plot2, ncol = 2,
top = textGrob("Fractional Brownian Motion Variance Estimation\n(1000 samples, n=100 each)",
gp = gpar(fontsize = 14, fontface = "bold")))
library(ggplot2)
# Function to generate Brownian bridge
generate_brownian_bridge <- function(n) {
# Standard Brownian motion
W <- c(0, cumsum(rnorm(n - 1, 0, sqrt(1/n))))
# Brownian bridge: B(t) = W(t) - t*W(1)
B <- W - seq(0, 1, length.out = n) * tail(W, 1)
return(B)
}
# Simulation parameters
n_trajectories <- 1000
n_points <- 100
time_points <- seq(0, 1, length.out = n_points)
# Generate trajectories
set.seed(123) # For reproducibility
trajectories <- matrix(NA, nrow = n_trajectories, ncol = n_points)
for (i in 1:n_trajectories) {
trajectories[i, ] <- generate_brownian_bridge(n_points)
}
# Calculate quantiles at each time point
lower_quantile <- apply(trajectories, 2, quantile, probs = 0.025)
upper_quantile <- apply(trajectories, 2, quantile, probs = 0.975)
# Create data frame for plotting
plot_data <- data.frame(
Time = rep(time_points, 3),
Value = c(colMeans(trajectories), lower_quantile, upper_quantile),
Type = rep(c("Mean", "2.5% Quantile", "97.5% Quantile"), each = n_points)
)
# Plot with ggplot2
ggplot(plot_data, aes(x = Time, y = Value, color = Type, linetype = Type)) +
geom_line(size = 0.8) +
scale_color_manual(values = c("black", "red", "red")) +
scale_linetype_manual(values = c("solid", "dotted", "dotted")) +
labs(title = "Brownian Bridge Trajectories with 95% Empirical Envelope",
subtitle = paste("Based on", n_trajectories, "simulated trajectories"),
x = "Time",
y = "Bridge Value",
color = "Statistic",
linetype = "Statistic") +
theme_minimal() +
theme(legend.position = "bottom",
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5),
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_blank()) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50", alpha = 0.5)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
library(ggplot2)
library(dplyr)
##
## Adjuntando el paquete: 'dplyr'
## The following object is masked from 'package:gridExtra':
##
## combine
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
# Function to generate Brownian bridge
generate_brownian_bridge <- function(n) {
W <- c(0, cumsum(rnorm(n - 1, 0, sqrt(1/n))))
B <- W - seq(0, 1, length.out = n) * tail(W, 1)
return(B)
}
# Simulation parameters
n_trajectories <- 1000
n_points <- 100
time_points <- seq(0, 1, length.out = n_points)
# Generate all trajectories
set.seed(123)
all_trajectories <- replicate(n_trajectories, generate_brownian_bridge(n_points))
# Prepare data for plotting (first 100 trajectories for visibility)
plot_data <- data.frame(
Time = rep(time_points, n_trajectories),
Value = as.vector(all_trajectories),
Trajectory = rep(1:n_trajectories, each = n_points)
) %>%
filter(Trajectory <= 100) # Plot only first 100 for clarity
# Calculate quantiles
quantile_data <- data.frame(
Time = time_points,
Lower = apply(all_trajectories, 1, quantile, probs = 0.025),
Upper = apply(all_trajectories, 1, quantile, probs = 0.975),
Mean = rowMeans(all_trajectories)
)
# Create the plot
ggplot() +
# Plot individual trajectories (light blue)
geom_line(data = plot_data,
aes(x = Time, y = Value, group = Trajectory),
color = "#ADD8E6", alpha = 0.3, size = 0.3) +
# Plot mean trajectory (dark blue)
geom_line(data = quantile_data,
aes(x = Time, y = Mean),
color = "blue", size = 0.8) +
# Plot quantiles (red dotted)
geom_line(data = quantile_data,
aes(x = Time, y = Lower),
color = "red", linetype = "dotted", size = 0.8) +
geom_line(data = quantile_data,
aes(x = Time, y = Upper),
color = "red", linetype = "dotted", size = 0.8) +
# Add horizontal reference line
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
# Labels and theme
labs(title = "Brownian Bridge with 95% Quantile Envelope",
subtitle = paste("Showing 100 of", n_trajectories, "simulated trajectories"),
x = "Time",
y = "Bridge Value",
caption = "Light blue: Individual trajectories\nDark blue: Mean trajectory\nRed dotted: 2.5% and 97.5% quantiles") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5),
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_blank(),
legend.position = "none")
