# Monte Carlo Simulation: Estimating Pi with Visualization
set.seed(42)
n_trials <- 5000
# 1. Generate random coordinates
x <- runif(n_trials, min = -1, max = 1)
y <- runif(n_trials, min = -1, max = 1)
# 2. Check if points are inside the unit circle (x^2 + y^2 <= 1)
inside_circle <- (x^2 + y^2) <= 1
pi_estimate <- 4 * sum(inside_circle) / n_trials
# 3. Calculate cumulative estimate to see convergence over time
cumulative_hits <- cumsum(inside_circle)
running_pi_est <- 4 * cumulative_hits / (1:n_trials)
# --- VISUALIZATION ---
# Set up a side-by-side plot layout
par(mfrow = c(1, 2))
# Plot 1: The Circle Distribution
plot(x, y, col = ifelse(inside_circle, "cornflowerblue", "firebrick1"),
pch = 20, cex = 0.5, asp = 1,
main = "Point Distribution",
sub = paste("Final Estimate:", round(pi_estimate, 4)))
# Draw the actual circle boundary for reference
symbols(0, 0, circles = 1, inches = FALSE, add = TRUE, lwd = 2)
# Plot 2: Convergence over Trials
plot(running_pi_est, type = "l", col = "darkgreen", lwd = 1.5,
main = "Convergence to Pi",
xlab = "Number of Trials", ylab = "Estimated Value")
abline(h = pi, col = "black", lty = 2, lwd = 2) # True value of Pi
legend("bottomright", legend = c("Estimate", "True Pi"),
col = c("darkgreen", "black"), lty = c(1, 2), bty = "n")
# Reset plotting parameters
par(mfrow = c(1, 1))
# Monte Carlo: Efficient Frontier Simulation
set.seed(123)
n_portfolios <- 10000
# 1. Input Assumptions (Expected Returns & Volatility)
# Asset A: High risk/return, B: Moderate, C: Low risk/return
exp_ret <- c(0.12, 0.08, 0.04)
volatility <- c(0.20, 0.12, 0.05)
cov_matrix <- diag(volatility^2) # Simplified: assuming no correlation
# 2. Simulation Loop
results <- data.frame(Risk = numeric(n_portfolios), Return = numeric(n_portfolios))
for (i in 1:n_portfolios) {
# Generate random weights that sum to 1
w <- runif(3)
w <- w / sum(w)
# Calculate Portfolio Return and Risk (Standard Deviation)
p_ret <- sum(w * exp_ret)
p_risk <- sqrt(t(w) %*% cov_matrix %*% w)
results[i, ] <- c(p_risk, p_ret)
}
# 3. Visualization
plot(results$Risk, results$Return, col = rgb(0.2, 0.4, 0.6, 0.2), pch = 20,
main = "Monte Carlo Efficient Frontier",
xlab = "Risk (Standard Deviation)", ylab = "Expected Return")
# Monte Carlo Simulation: Project Schedule Risk
set.seed(456)
n_simulations <- 10000
# 1. Define task durations based on different probability distributions
# Task A: Normal distribution
task_a <- rnorm(n_simulations, mean = 10, sd = 2)
# Task B: Triangular distribution (using a simple logic for min, mode, max)
# Formula: runif() transformed or using the 'triangle' package logic
# Here we'll simulate a Build phase with a bias toward 10 days
task_b <- replicate(n_simulations, {
u <- runif(1)
if(u < 0.33) runif(1, 5, 10) else runif(1, 10, 20)
})
# Task C: Uniform distribution
task_c <- runif(n_simulations, min = 8, max = 12)
# 2. Calculate Total Project Duration
total_duration <- task_a + task_b + task_c
# 3. Calculate the Probability of finishing within 35 days
prob_on_time <- mean(total_duration <= 35)
# --- VISUALIZATION ---
# Set up plotting area
par(mfrow = c(1, 1))
# Histogram of results
hist(total_duration, breaks = 50, col = "skyblue", border = "white",
main = "Distribution of Total Project Duration",
xlab = "Total Days", ylab = "Frequency")
# Add a vertical line for the 35-day deadline
abline(v = 35, col = "red", lwd = 3, lty = 2)
# Add text results to the plot
text(42, 400, paste("Prob <= 35 days:", round(prob_on_time * 100, 1), "%"),
col = "red", font = 2)
# Monte Carlo: Full Portfolio Analytics & Frontier
set.seed(42)
n_sims <- 10000
# 1. SETUP: Inputs (Tech, Energy, Gold, Bonds)
assets <- c("Tech", "Energy", "Gold", "Bonds")
exp_ret <- c(0.14, 0.10, 0.04, 0.02)
volatility <- c(0.25, 0.20, 0.15, 0.05)
risk_free <- 0.03
# Covariance Matrix construction
cor_mat <- matrix(c(1.0, 0.4, -0.1, 0.1, 0.4, 1.0, 0.2, 0.0, -0.1, 0.2, 1.0, 0.3, 0.1, 0.0, 0.3, 1.0), nrow=4)
cov_mat <- diag(volatility) %*% cor_mat %*% diag(volatility)
# 2. MONTE CARLO SIMULATION
results <- data.frame(Return = numeric(n_sims), Risk = numeric(n_sims), Sharpe = numeric(n_sims))
weights_matrix <- matrix(nrow = n_sims, ncol = length(assets))
for (i in 1:n_sims) {
w <- runif(4); w <- w / sum(w)
p_ret <- sum(w * exp_ret)
p_risk <- sqrt(t(w) %*% cov_mat %*% w)
p_sharpe <- (p_ret - risk_free) / p_risk
results[i, ] <- c(p_ret, p_risk, p_sharpe)
weights_matrix[i, ] <- w
}
# 3. IDENTIFY KEY PORTFOLIOS
max_sharpe_idx <- which.max(results$Sharpe)
min_risk_idx <- which.min(results$Risk)
# 4. PRINTING CALCULATED VALUES
cat("--- TANGENCY PORTFOLIO (MAX SHARPE) ---\n")
## --- TANGENCY PORTFOLIO (MAX SHARPE) ---
cat("Portfolio Return: ", round(results$Return[max_sharpe_idx] * 100, 2), "%\n")
## Portfolio Return: 10.81 %
cat("Portfolio Risk (Volatility): ", round(results$Risk[max_sharpe_idx] * 100, 2), "%\n")
## Portfolio Risk (Volatility): 16.17 %
cat("Sharpe Ratio: ", round(results$Sharpe[max_sharpe_idx], 4), "\n\n")
## Sharpe Ratio: 0.4828
cat("--- MINIMUM VARIANCE PORTFOLIO ---\n")
## --- MINIMUM VARIANCE PORTFOLIO ---
cat("Portfolio Return: ", round(results$Return[min_risk_idx] * 100, 2), "%\n")
## Portfolio Return: 2.55 %
cat("Portfolio Risk: ", round(results$Risk[min_risk_idx] * 100, 2), "%\n\n")
## Portfolio Risk: 4.99 %
# --- VISUALIZATION ---
par(mfrow = c(1, 2), mar = c(7, 4, 4, 2)) # Increased bottom margin for labels
# Graph 1: Efficient Frontier & Capital Market Line (CML)
plot(results$Risk, results$Return, col = rgb(0.1, 0.2, 0.5, 0.1), pch = 20,
main = "Efficient Frontier & CML", xlab = "Risk (Std Dev)", ylab = "Return")
abline(a = risk_free, b = results$Sharpe[max_sharpe_idx], col = "darkgreen", lwd = 2)
points(results$Risk[max_sharpe_idx], results$Return[max_sharpe_idx], col = "red", pch = 19, cex = 1.5)
points(results$Risk[min_risk_idx], results$Return[min_risk_idx], col = "gold", pch = 19, cex = 1.5)
# Graph 2: Asset Allocation (Max Sharpe Portfolio)
# las=2 rotates labels; cex.names shrinks them if needed
bp <- barplot(weights_matrix[max_sharpe_idx, ], names.arg = assets, col = "steelblue",
main = "Optimal Asset Allocation", ylab = "Weight %", ylim = c(0, 1.1),
las = 2, cex.names = 0.8)
# Add percentage labels above bars
text(x = bp, y = weights_matrix[max_sharpe_idx, ] + 0.05,
labels = paste0(round(weights_matrix[max_sharpe_idx, ] * 100, 1), "%"), cex = 0.8)
