library(tidyverse)
library(tidyquant)
library(patchwork)
library(scales)
<- c("PG", "JPM", "NVDA")
symbols <- tq_get(symbols, get = "stock.prices", from = "2012-12-31", to = "2022-12-31")
mult_stocks $symbol <- mult_stocks$symbol <- factor(mult_stocks$symbol, levels = c("PG", "JPM", "NVDA"))
mult_stocks
# tq_mutate_fun_options() returns list of compatible mutate functions by pkg
<- mult_stocks |>
all_returns_daily group_by(symbol) |>
tq_transmute(select = adjusted,
mutate_fun = periodReturn,
period = "daily",
type = "log")
<- mult_stocks |>
all_returns_monthly group_by(symbol) |>
tq_transmute(select = adjusted,
mutate_fun = periodReturn,
period = "monthly",
type = "log")
# reframe() has apparently replaced summarize()
<- all_returns_daily |>
quantiles_daily group_by(symbol) |>
reframe(quantiles = quantile(daily.returns, probs = c(0.05, 0.95)))
<- all_returns_monthly |>
quantiles_monthly group_by(symbol) |>
reframe(quantiles = quantile(monthly.returns, probs = c(0.05, 0.95)))
# 5% quantile for each stock, DAILY
<- quantiles_daily$quantiles[1]
PG_05d <- quantiles_daily$quantiles[3]
JPM_05d <- quantiles_daily$quantiles[5]
NVDA_05d <- mean(all_returns_daily$daily.returns)
mean_d
# 5% quantile for each stock, MONTHLY
<- quantiles_monthly$quantiles[1]
PG_05m <- quantiles_monthly$quantiles[3]
JPM_05m <- quantiles_monthly$quantiles[5]
NVDA_05m <- mean(all_returns_monthly$monthly.returns)
mean_m
# I probably spend too much time tinkering with colors
# col_ticker_fills <- c("PG" = "blue", "JPM" = "yellow", "NVDA" = "red")
# col_ticker_fills <- c("PG" = "#90EE90", "JPM" = "#ff6347", "NVDA" = "#8B0000")
<- c("PG" = "chartreuse2", "JPM" = "dodgerblue2", "NVDA" = "coral1")
col_ticker_fills <- c("PG" = "chartreuse3", "JPM" = "dodgerblue3", "NVDA" = "coral3")
col_ticker_colors <- "chartreuse3"; col_JPM_line <- "dodgerblue3"; col_NVDA_line <- "coral3"
col_PG_line
<- all_returns_daily |>
p_hist_daily ggplot(aes(x = daily.returns, fill = symbol, color = symbol)) +
geom_density(alpha = 0.50) +
geom_vline(xintercept = PG_05d, color = col_PG_line, linetype = "dashed", linewidth = 1) +
geom_vline(xintercept = JPM_05d, color = col_JPM_line, linetype = "dashed", linewidth = 1) +
geom_vline(xintercept = NVDA_05d, color = col_NVDA_line, linetype = "dashed", linewidth = 1) +
geom_vline(xintercept = mean_d, color = "black", linewidth = 0.4) +
scale_fill_manual(values = col_ticker_fills) +
scale_color_manual(values = col_ticker_fills) +
theme_minimal() +
coord_cartesian(xlim = c(-0.25, 0.25)) +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
<- all_returns_monthly |>
p_hist_monthly ggplot(aes(x = monthly.returns, fill = symbol, color = symbol)) +
geom_density(alpha = 0.50) +
geom_vline(xintercept = PG_05m, color = col_PG_line, linetype = "dashed", linewidth = 1) +
geom_vline(xintercept = JPM_05m, color = col_JPM_line, linetype = "dashed", linewidth = 1) +
geom_vline(xintercept = NVDA_05m, color = col_NVDA_line, linetype = "dashed", linewidth = 1) +
geom_vline(xintercept = mean_m, color = "black", linewidth = 0.4) +
scale_fill_manual(values = col_ticker_fills) +
scale_color_manual(values = col_ticker_fills) +
theme_minimal() +
coord_cartesian(xlim = c(-0.25, 0.25)) +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
<- p_hist_daily / p_hist_monthly
p_hs p_hs
Contents
- Historical Simulation: Basic and Bootstrap
- Monte Carlo
- Parametric; aka, analytical
Historical simulation (HS)
Basic HS
Bootstrap HS
First let’s simulate a single one-year forward (+ 12 months) path
# all_returns_monthly is tidy, but the simulation wants a wide format
<- all_returns_monthly |>
all_returns_monthly_wide pivot_wider(names_from = symbol, values_from = monthly.returns)
# Initial investment (into each stock)
<- 100
initial_investment <- setNames(data.frame(t(rep(initial_investment, times = length(symbols)))), symbols)
portfolio <- 12 # Set X months for the simulation
months_to_simulate
# Simulate one forward month
<- function(portfolio, historical_returns_wide) {
simulate_one_month # Randomly sample one month's returns (with replacement)
<- historical_returns_wide |>
sampled_returns sample_n(1, replace = TRUE) |> select(-date)
# Apply the sampled log returns to the current portfolio value
<- portfolio * exp(sampled_returns)
updated_portfolio # FOR TESTING: print(sampled_returns[1,]); print(updated_portfolio[1,])
return(updated_portfolio)
}
# Run the simulation for X months
set.seed(123) # For reproducibility of results
<- tibble(Month = 0, TotalValue = sum(portfolio))
simulation_results for (i in 1:months_to_simulate) {
<- simulate_one_month(portfolio, all_returns_monthly_wide)
portfolio <- simulation_results |>
simulation_results add_row(Month = i, TotalValue = sum(portfolio))
}
print(simulation_results)
# A tibble: 13 × 2
Month TotalValue
<dbl> <dbl>
1 0 300
2 1 294.
3 2 325.
4 3 329.
5 4 316.
6 5 314.
7 6 349.
8 7 356.
9 8 312.
10 9 313.
11 10 326.
12 11 349.
13 12 336.
Next let’s add a loop to run the simulation multiple times; e.g. 20 trials
library(RColorBrewer)
# library(ggside) couldn't pull it off!
# Function to simulate portfolio over X months, where each trail has a trial_id
<- function(months_to_simulate, historical_returns_wide, initial_investment, trial_id) {
simulate_portfolio # Initialize portfolio
<- setNames(data.frame(t(rep(initial_investment, times = length(symbols)))), symbols)
portfolio # Initialize results data frame with Month 0
<- tibble(Month = 0, TotalValue = sum(portfolio), Trial = as.factor(trial_id))
simulation_results for (i in 1:months_to_simulate) {
<- simulate_one_month(portfolio, historical_returns_wide)
portfolio <- simulation_results |>
simulation_results add_row(Month = i, TotalValue = sum(portfolio), Trial = as.factor(trial_id))
}return(simulation_results)
}
<- 12
months_to_simulate <- 20 # Number of trials
num_trials
set.seed(123)
<- map_df(1:num_trials,
all_trials ~simulate_portfolio(months_to_simulate, all_returns_monthly_wide, initial_investment, .x),
.id = "Trial_ID")
<- all_trials |>
final_month_values_df filter(Month == max(all_trials$Month))
# Plot the results using ggplot2
<- ggplot(all_trials, aes(x = Month, y = TotalValue, group = Trial, color = Trial)) +
p_forward_sim geom_line() +
scale_color_viridis_d(option = "plasma", direction = -1) +
theme_minimal() +
scale_x_continuous(breaks = 1:12, limits = c(0,12)) +
labs(x = "Month",
y = "Portfolio Value") +
theme(legend.position = "none")
<- final_month_values_df |> ggplot(aes(x = TotalValue)) +
density_plot geom_density(fill = "#933fbd", alpha = 0.5) +
theme_minimal() +
theme(
axis.title = element_blank(),
axis.text = element_blank()
+
) coord_flip()
<- p_forward_sim + density_plot +
p_boot plot_layout(ncol = 2, widths = c(3, 1))
p_boot
<- final_month_values_df |>
final_month_values_vct pull(TotalValue) # 'pull' extracts the column as a vector
# Calculate and print the quantiles for the final month's values
<- quantile(final_month_values_vct, probs = c(0, 0.01, 0.05, 0.50, 1.0))
quantiles_final_month print(quantiles_final_month)
0% 1% 5% 50% 100%
218.1581 233.4910 294.8226 371.7791 567.0854
Monte carlos simulation (MCS)
First a single path:
<- function(price, mu, sigma, periods) {
sim_gbm <- 1/52 # we'll assume period is one week
dt <- numeric(periods)
prices 1] <- price
prices[
for(t in 2:periods) {
<- rnorm(1) # random standar normal quantile
Z
# GBM
<- prices[t-1]*exp((mu -0.5*sigma^2)*dt +sigma*sqrt(dt)*Z)
prices[t]
}return(prices)
}
set.seed(952347)
<- sim_gbm(100, 0.09, 0.3, 50)
sim_prices <- as_tibble(sim_prices) |> rownames_to_column("period")
sim_prices_df $period <- as.integer(sim_prices_df$period)
sim_prices_df
# Plot the simulated stock price path
|>
sim_prices_df ggplot(aes(x = period, y = value)) +
geom_line() +
theme_minimal()
Now multiple (say 20) trials:
<- function(price, mu, sigma, periods, simulations) {
sim_gbm_matrix <- 1/52
dt <- matrix(price, nrow = periods, ncol = simulations)
prices for(t in 2:periods) {
<- rnorm(simulations)
Z <- prices[t-1, ] * exp((mu - 0.5 * sigma^2) * dt + sigma * sqrt(dt) * Z)
prices[t, ]
}return(prices)
}
set.seed(84923)
<- 20
simulations <- sim_gbm_matrix(100, 0.09, 0.3, 50, simulations)
sim_prices_matrix
<- as_tibble(sim_prices_matrix, .name_repair = "minimal")
sim_prices_df names(sim_prices_df) <- 1:ncol(sim_prices_df)
<- sim_prices_df |> rownames_to_column("period")
sim_prices_df $period <- as.integer(sim_prices_df$period)
sim_prices_df
<- sim_prices_df |>
sim_prices_long pivot_longer(cols = !period, names_to = "Trial", values_to = "Price")
<- sim_prices_long |> ggplot(aes(x = period, y = Price, group = Trial, color = Trial)) +
p_mcs geom_line() +
labs(x = "Week",
y = "Price") +
theme_minimal() +
theme(
legend.position = "none"
+
) scale_color_viridis_d()
+ ggtitle("GBM MCS: 20 trials x 50 weeks") p_mcs
Parametric: normally distributed arithmetic returns
The basic starting point: normal, arithmetic returns and scaling per the square root rule (SRR) which asumes i.i.d. returns.
library(ggplot2)
# Parameters
<- 0.09
mu_pa <- 0.20
sigma_pa <- 12/52
horizon
<- mu_pa * horizon # horizon return
mu <- sigma_pa * sqrt(horizon) # horizon volatility
sigma <- 0.95
confidence_level <- 100 # Initial investment
investment_value
# Find the Z-score that corresponds to the confidence level
<- qnorm(confidence_level)
z_score <- (-mu + sigma * z_score) * investment_value
VaR
# Create a dataframe of returns for plotting
<- data.frame(Returns = seq(-0.30, 0.30, by = 0.001))
returns $Density <- dnorm(returns$Returns, mean = mu, sd = sigma)
returns
# Plot with ggplot
<- ggplot(returns, aes(x = Returns, y = Density)) +
p_pVaR geom_line(color = "dodgerblue", linewidth = 1.3) +
geom_vline(xintercept = -VaR / investment_value, color = "darkred", linetype = "dashed", linewidth = 1.3) +
geom_area(data = subset(returns, Returns < -VaR / investment_value),
aes(x = Returns, y = Density), fill = "coral1", alpha = 0.5) +
geom_vline(xintercept = mu, color = "chartreuse3", linetype = "dashed", linewidth = 1.3) +
geom_vline(xintercept = 0, color = "black", linewidth = 0.4) +
scale_x_continuous(labels = percent_format(),
breaks = seq(-0.30, 0.30, by = 0.1)) +
theme_minimal() +
theme(
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank()
+
) annotate("text", x = -VaR / investment_value + .02,
y = max(returns$Density)/3,
label = paste("VaR at", scales::percent(confidence_level), "=", round(VaR, 3)),
hjust = 1.2, vjust = 0, size = 4, color = "darkred", fontface = "bold")
+ ggtitle("Normal, arithmetic i.i.d. returns") p_pVaR
For fun: Using patchwork to print the social thumbnail!
# layout <- "
# AABBCC
# AABBCC
# DDDDEE
# DDDDEE
# "
#
# p_hs + p_pVaR + p_boot + p_mcs +
# plot_layout(design = layout)
/ p_boot / (p_mcs + p_pVaR) p_hs