# Load packages

# Core
library(tidyverse)
library(tidyquant)

# Source function
source("../00_scripts/simulate_accumulation.R")

1 Import stock prices

Revise the code below.

symbols <- c("RTX", "GD", "LMT", "BA")

prices <- tq_get(x    = symbols,
                 get  = "stock.prices",    
                 from = "2012-12-31",
                 to   = "2022-11-30")

2 Convert prices to returns

asset_returns_tbl <- prices %>%
    
    group_by(symbol) %>%
    
    tq_transmute(select     = adjusted, 
                 mutate_fun = periodReturn, 
                 period     = "monthly",
                 type       = "log") %>%
    
    slice(-1) %>%
    
    ungroup() %>%
    
    set_names(c("asset", "date", "returns"))

3 Assign a weight to each asset

Revise the code for weights.

# symbols
symbols <- asset_returns_tbl %>% distinct(asset) %>% pull()
symbols
## [1] "BA"  "GD"  "LMT" "RTX"
# weights
weights <- c(0.35, 0.30, 0.20, 0.15)
weights
## [1] 0.35 0.30 0.20 0.15
w_tbl <- tibble(symbols, weights)
w_tbl
## # A tibble: 4 × 2
##   symbols weights
##   <chr>     <dbl>
## 1 BA         0.35
## 2 GD         0.3 
## 3 LMT        0.2 
## 4 RTX        0.15

4 Build a portfolio

portfolio_returns_tbl <- asset_returns_tbl %>%
    
    tq_portfolio(assets_col = asset, 
                 returns_col = returns, 
                 weights = w_tbl, 
                 rebalance_on = "months", 
                 col_rename = "returns")

portfolio_returns_tbl
## # A tibble: 119 × 2
##    date       returns
##    <date>       <dbl>
##  1 2013-01-31 -0.0224
##  2 2013-02-28  0.0349
##  3 2013-03-28  0.0727
##  4 2013-04-30  0.0405
##  5 2013-05-31  0.0642
##  6 2013-06-28  0.0184
##  7 2013-07-31  0.0764
##  8 2013-08-30 -0.0114
##  9 2013-09-30  0.0773
## 10 2013-10-31  0.0423
## # … with 109 more rows

5 Simulating growth of a dollar

# Get mean portfolio return
mean_port_return <- mean(portfolio_returns_tbl$returns)
mean_port_return
## [1] 0.01116809
# Get standard deviation of portfolio returns
stddev_port_return <- sd(portfolio_returns_tbl$returns)
stddev_port_return
## [1] 0.06662271

6 Simulation function

No need

7 Running multiple simulations

# Create a vector of 1s as a starting point
sims <- 51
starts <- rep(100, sims) %>% 
    set_names(paste0("sims", 1:sims))

starts
##  sims1  sims2  sims3  sims4  sims5  sims6  sims7  sims8  sims9 sims10 sims11 
##    100    100    100    100    100    100    100    100    100    100    100 
## sims12 sims13 sims14 sims15 sims16 sims17 sims18 sims19 sims20 sims21 sims22 
##    100    100    100    100    100    100    100    100    100    100    100 
## sims23 sims24 sims25 sims26 sims27 sims28 sims29 sims30 sims31 sims32 sims33 
##    100    100    100    100    100    100    100    100    100    100    100 
## sims34 sims35 sims36 sims37 sims38 sims39 sims40 sims41 sims42 sims43 sims44 
##    100    100    100    100    100    100    100    100    100    100    100 
## sims45 sims46 sims47 sims48 sims49 sims50 sims51 
##    100    100    100    100    100    100    100
# Simulate
# For reproducible research
set.seed(1234)

monte_carlo_sim_51 <- starts %>% 
    
    # Simulate
    map_dfc(.x = ., 
            .f = ~simulate_accumulation(initial_value = .x, 
                                        N             = 240, 
                                        mean_return   = mean_port_return, 
                                        sd_return     = stddev_port_return)) %>% 
    
    # Add column month
    mutate(month = 1:nrow(.)) %>% 
    select(month, everything()) %>% 
    
    # Rearrange column names
    set_names(c("month", names(starts))) %>% 
    
    # Transform to long form
    pivot_longer(cols = -month, names_to = "sim", values_to = "growth")

monte_carlo_sim_51
## # A tibble: 12,291 × 3
##    month sim    growth
##    <int> <chr>   <dbl>
##  1     1 sims1     100
##  2     1 sims2     100
##  3     1 sims3     100
##  4     1 sims4     100
##  5     1 sims5     100
##  6     1 sims6     100
##  7     1 sims7     100
##  8     1 sims8     100
##  9     1 sims9     100
## 10     1 sims10    100
## # … with 12,281 more rows
# Find quantiles
monte_carlo_sim_51 %>% 
    
    group_by(sim) %>% 
    summarise(growth = last(growth)) %>% 
    ungroup() %>%
    pull(growth) %>% 
    
    quantile(probs = c(0, 0.25, 0.5, 0.75, 1)) %>% 
    round(2)
##      0%     25%     50%     75%    100% 
##   99.92  566.07 1165.56 2125.22 5339.70

8 Visualizing simulations with ggplot

Line Plot of Simulations with Max, Median, and Min

# Step 1: Summarize data into max, median, and min of the last value
sim_summary <- monte_carlo_sim_51 %>% 
    
    group_by(sim) %>% 
    summarise(growth = last(growth)) %>% 
    ungroup() %>% 
    
    summarise(max    = max(growth), 
              median = median(growth), 
              min    = min(growth)) 

sim_summary
## # A tibble: 1 × 3
##     max median   min
##   <dbl>  <dbl> <dbl>
## 1 5340.  1166.  99.9
# Step 2 Plot
monte_carlo_sim_51 %>% 
    
    # Filter for max, median, and min sim
    group_by(sim) %>% 
    filter(last(growth) == sim_summary$max | 
               last(growth) == sim_summary$median | 
               last(growth) == sim_summary$min) %>% 
    ungroup() %>% 
    
    # Plot
    ggplot(aes(x = month, y = growth, color = sim)) + 
    geom_line() + 
    
    theme(legend.position = "none") + 
    theme(plot.title = element_text(hjust = 0.5)) + 
    theme(plot.subtitle = element_text(hjust = 0.5)) + 
    
    labs(title = "Simulating growth of $100 over 120 months", 
         subtitle = "Maximum, Median, and Minimum Simulation")

Based on the Monte Carlo simulation results, how much should you expect from your $100 investment after 20 years? What is the best-case scenario? What is the worst-case scenario? What are limitations of this simulation analysis?

I would expect around the median amount of $1,166 after 20 years. The best case scenario of $5,340 seems a bit high but it could be possible. The worst-case scenario say that this portfolio would lose 1 cent over 20 years. This could happen depending on world events. The limitations of this simulation is that there is no way of actually knowing what could happen. The furure is unknown. This simulation is good for approximating what could happen but I would probably want to run more than 51 simulations. It is really intersting to be able to run a simulation like this. It would absolutely help when making decisions for portfolios. And it is easy to adjust now that it is set up.