# 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("AAPL", "MMM", "MSFT", "AMZN", "TSLA")

prices <- tq_get(x    = symbols,
                 get  = "stock.prices",    
                 from = "2012-12-31",
                 to   = "2023-04-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] "AAPL" "AMZN" "MMM"  "MSFT" "TSLA"
# weights
weights <- c(0.25, 0.25, 0.2, 0.2, 0.1)
weights
## [1] 0.25 0.25 0.20 0.20 0.10
w_tbl <- tibble(symbols, weights)
w_tbl
## # A tibble: 5 × 2
##   symbols weights
##   <chr>     <dbl>
## 1 AAPL       0.25
## 2 AMZN       0.25
## 3 MMM        0.2 
## 4 MSFT       0.2 
## 5 TSLA       0.1

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: 124 × 2
##    date        returns
##    <date>        <dbl>
##  1 2013-01-31  0.00688
##  2 2013-02-28 -0.00282
##  3 2013-03-28  0.0214 
##  4 2013-04-30  0.0494 
##  5 2013-05-31  0.103  
##  6 2013-06-28 -0.0181 
##  7 2013-07-31  0.0737 
##  8 2013-08-30  0.0311 
##  9 2013-09-30  0.0440 
## 10 2013-10-31  0.0650 
## # … with 114 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.01743237
# Get standard deviation of portfolio returns
stddev_port_return <- sd(portfolio_returns_tbl$returns)
stddev_port_return
## [1] 0.06093678

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(paste("sim", 1:sims, sep = ""))
starts
##  sim1  sim2  sim3  sim4  sim5  sim6  sim7  sim8  sim9 sim10 sim11 sim12 sim13 
##   100   100   100   100   100   100   100   100   100   100   100   100   100 
## sim14 sim15 sim16 sim17 sim18 sim19 sim20 sim21 sim22 sim23 sim24 sim25 sim26 
##   100   100   100   100   100   100   100   100   100   100   100   100   100 
## sim27 sim28 sim29 sim30 sim31 sim32 sim33 sim34 sim35 sim36 sim37 sim38 sim39 
##   100   100   100   100   100   100   100   100   100   100   100   100   100 
## sim40 sim41 sim42 sim43 sim44 sim45 sim46 sim47 sim48 sim49 sim50 sim51 
##   100   100   100   100   100   100   100   100   100   100   100   100
# Simulate
monte_carlo_sim_51 <- starts %>%
    # Simulate
    map_dfc(simulate_accumulation,
            N     = 240,
            mean  = mean_port_return,
            stdev = stddev_port_return) %>%
    # Add the column, month
    mutate(month = seq(1:nrow(.))) %>%
    # Arrange column names
    select(month, everything()) %>%
    set_names(c("month", names(starts))) %>%
    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 sim1     100
##  2     1 sim2     100
##  3     1 sim3     100
##  4     1 sim4     100
##  5     1 sim5     100
##  6     1 sim6     100
##  7     1 sim7     100
##  8     1 sim8     100
##  9     1 sim9     100
## 10     1 sim10    100
## # … with 12,281 more rows
# Calculate the quantiles for simulated values
probs <- c(.005, .025, .25, .5, .75, .975, .995)
monte_carlo_sim_51 %>%
    group_by(sim) %>%
    summarise(growth = last(growth)) %>%
    ungroup() %>%
    pull(growth) %>%
    # Find the quantiles
    quantile(probs = probs) %>%
    round(2)
##     0.5%     2.5%      25%      50%      75%    97.5%    99.5% 
##   843.58  1379.15  3079.90  4941.04  9390.89 21136.37 43503.04

8 Visualizing simulations with ggplot

monte_carlo_sim_51 %>%
    ggplot(aes(x = month, y = growth, col = sim)) +
    geom_line() +
    theme(legend.position = "none")

# Simplify the plot
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 50864.  4941.  691.
monte_carlo_sim_51 %>%
    group_by(sim) %>%
    filter(last(growth) == sim_summary$max |
           last(growth) == sim_summary$median |
           last(growth) == sim_summary$min) %>%
    # Plot
    ggplot(aes(month, growth, col = sim)) +
    geom_line() +
    theme()

Line Plot of Simulations with Max, Median, and Min

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?

My 100 dollar investment has a very wide range of expected returns but the best case scenario it will be worth over $10,000 and the worst case scenario is to small to see since the best case is so high. Some limitations of this analysis is that my companies have experienced very high growth in the past and that growth is not possible to keep going for the future.