# Load packages

# Core
library(tidyverse)
library(tidyquant)

1 Import stock prices

Revise the code below.

symbols <- c("NVDA", "TSLA", "AMD", "MSFT")

prices <- tq_get(x    = symbols,
                 get  = "stock.prices",    
                 from = "2012-12-31",
                 to   = "2024-07-01")

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] "AMD"  "MSFT" "NVDA" "TSLA"
# weights
weights <- c(0.25, 0.25, 0.25, 0.25)
weights
## [1] 0.25 0.25 0.25 0.25
w_tbl <- tibble(symbols, weights)
w_tbl
## # A tibble: 4 × 2
##   symbols weights
##   <chr>     <dbl>
## 1 AMD        0.25
## 2 MSFT       0.25
## 3 NVDA       0.25
## 4 TSLA       0.25

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: 138 × 2
##    date       returns
##    <date>       <dbl>
##  1 2013-01-31  0.0524
##  2 2013-02-28 -0.0146
##  3 2013-03-28  0.0375
##  4 2013-04-30  0.168 
##  5 2013-05-31  0.264 
##  6 2013-06-28  0.0182
##  7 2013-07-31  0.0229
##  8 2013-08-30  0.0422
##  9 2013-09-30  0.0844
## 10 2013-10-31 -0.0709
## # ℹ 128 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.03222546
# Get standard deviation of portfolio returns
stddev_port_return <- sd(portfolio_returns_tbl$returns)
stddev_port_return
## [1] 0.09495641
# Construct a normal distribution
simulated_monthly_returns <- rnorm(240, mean_port_return, stddev_port_return)
simulated_monthly_returns
##   [1]  0.094053650  0.019808055  0.125149478 -0.008840974  0.088643187
##   [6] -0.046105613  0.012606941  0.131319064  0.125302020  0.084256861
##  [11] -0.063968451  0.021213585  0.045848342  0.129821337  0.149062600
##  [16]  0.036295972  0.067768308 -0.045872087  0.059016956  0.031462485
##  [21] -0.011595529  0.143994016  0.189974799 -0.048579774 -0.118744873
##  [26]  0.108338723 -0.007965295  0.024983648  0.034207521  0.067978837
##  [31]  0.048872750  0.177775949  0.062233331 -0.193151711  0.067841443
##  [36]  0.150367418 -0.122707625  0.128323842  0.050705444 -0.042185559
##  [41]  0.124164004  0.150939072  0.002567559  0.013780392  0.341273962
##  [46]  0.096121173 -0.214490311  0.153715324  0.216114246 -0.068709767
##  [51]  0.059268153 -0.065287260  0.185190550  0.013794194 -0.143964969
##  [56] -0.287643946 -0.019583085 -0.060938229 -0.036304017  0.091833744
##  [61]  0.006000159  0.096594870  0.085666436  0.068370320 -0.122733748
##  [66]  0.047869595 -0.034500768  0.140058046  0.056452401 -0.019636182
##  [71]  0.036418869  0.149565300  0.018644002 -0.051780386 -0.027364773
##  [76]  0.088587179 -0.055198576 -0.030173526  0.171176827  0.175343873
##  [81]  0.076927239 -0.026144344  0.244472701  0.147819791  0.031021314
##  [86]  0.214080068 -0.122164851  0.104826536  0.087514387  0.174809037
##  [91] -0.007168830 -0.142701922  0.064683564 -0.026111931  0.140076350
##  [96]  0.061790784  0.082277238  0.199422399 -0.008747502  0.117976955
## [101]  0.121128197 -0.008535493  0.060617731 -0.080214824  0.036823900
## [106]  0.119367173  0.100205713  0.030908185  0.116962701  0.020096244
## [111]  0.081532090  0.120442440  0.061450401 -0.065205644  0.034941075
## [116] -0.160832700 -0.046364409 -0.017639601  0.134819587 -0.177192342
## [121] -0.050525512  0.144179479 -0.010935734  0.091073344  0.030759614
## [126]  0.062302780 -0.036785521  0.015570448  0.034491074  0.037198328
## [131]  0.100712039 -0.032288991  0.162598061  0.215655648  0.011498514
## [136]  0.185168490 -0.107355281  0.090473242 -0.117689526 -0.226494172
## [141]  0.036313321 -0.051490883  0.066249013 -0.039491514  0.072945953
## [146]  0.077039986  0.035261854  0.041473369 -0.122832257 -0.039111399
## [151] -0.029629194 -0.006119228  0.031263186 -0.086466656  0.004019295
## [156] -0.086422539  0.107933010  0.032680855  0.186655238 -0.033598526
## [161]  0.014860498  0.197967804 -0.044334759  0.005050655 -0.054178283
## [166]  0.021623129  0.084928373  0.003347401  0.213715942  0.061199520
## [171]  0.037237945 -0.070962843 -0.171545315 -0.054837835 -0.020215344
## [176]  0.060209077  0.010081407  0.210649584 -0.060342937  0.023696685
## [181]  0.204954287 -0.003001609 -0.060783308  0.126888524 -0.194271971
## [186] -0.061569308  0.187092603  0.074253031  0.228772868  0.017093128
## [191]  0.056761368  0.037749175 -0.015797485  0.109514270  0.089280015
## [196]  0.175636862  0.054859896  0.029442395  0.077199733  0.030127447
## [201]  0.239725705 -0.010033773  0.064652309 -0.141958180 -0.061068135
## [206]  0.029693746  0.009066125  0.127257126  0.182519193 -0.098852419
## [211] -0.095567560  0.062389005 -0.010883762  0.130230000  0.094636782
## [216] -0.037293727 -0.006250113 -0.060816572  0.067859357  0.011568523
## [221] -0.062032254 -0.083290244  0.166884018  0.039398822  0.094160680
## [226] -0.002673206  0.102970354  0.018427358  0.120900093 -0.078376788
## [231]  0.118884966  0.180505254 -0.023630253  0.057923463  0.061440730
## [236] -0.118463860  0.164292953  0.031448642 -0.082237436  0.094430058
# Add a dollar
simulated_returns_add_1 <- tibble(returns = c(1, 1 + simulated_monthly_returns))
simulated_returns_add_1
## # A tibble: 241 × 1
##    returns
##      <dbl>
##  1   1    
##  2   1.09 
##  3   1.02 
##  4   1.13 
##  5   0.991
##  6   1.09 
##  7   0.954
##  8   1.01 
##  9   1.13 
## 10   1.13 
## # ℹ 231 more rows
# Calculate the cumulative growth of a dollar
simulated_growth <- simulated_returns_add_1 %>%
    mutate(growth = accumulate(returns, function(x, y) x*y)) %>%
    select(growth)

simulated_growth
## # A tibble: 241 × 1
##    growth
##     <dbl>
##  1   1   
##  2   1.09
##  3   1.12
##  4   1.26
##  5   1.24
##  6   1.35
##  7   1.29
##  8   1.31
##  9   1.48
## 10   1.67
## # ℹ 231 more rows
# Check the compound annual growth rate
cagr <- ((simulated_growth$growth[nrow(simulated_growth)]^(1/10)) - 1) * 100
cagr
## [1] 111.1729

6 Simulation function

simulate_accmulation <- function(initial_value, N, mean_return, sd_return) {                
    
    # Add a dollar
    simulated_returns_add_1 <- tibble(returns = c(initial_value, 1 + rnorm(N, mean_return, sd_return)))
   
    # Calculate the cumulative growth of a dollar
    simulated_growth <- simulated_returns_add_1 %>%
        mutate(growth = accumulate(returns, function(x, y) x*y)) %>%
        select(growth)

    return(simulated_growth)
}
    
simulate_accmulation(initial_value = 100, N = 240, mean_return = 0.005, sd_return = 0.01) %>%
    tail()
## # A tibble: 6 × 1
##   growth
##    <dbl>
## 1   339.
## 2   343.
## 3   342.
## 4   345.
## 5   353.
## 6   358.

7 Running multiple simulations

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

starts
##  sim1  sim2  sim3  sim4  sim5  sim6  sim7  sim8  sim9 sim10 sim11 sim12 sim13 
##     1     1     1     1     1     1     1     1     1     1     1     1     1 
## sim14 sim15 sim16 sim17 sim18 sim19 sim20 sim21 sim22 sim23 sim24 sim25 sim26 
##     1     1     1     1     1     1     1     1     1     1     1     1     1 
## sim27 sim28 sim29 sim30 sim31 sim32 sim33 sim34 sim35 sim36 sim37 sim38 sim39 
##     1     1     1     1     1     1     1     1     1     1     1     1     1 
## sim40 sim41 sim42 sim43 sim44 sim45 sim46 sim47 sim48 sim49 sim50 sim51 
##     1     1     1     1     1     1     1     1     1     1     1     1
# Simulate 
Portfolio_sim <- starts %>%
    
    # Simulate
    map_dfc(.x = .,
            .f = ~simulate_accmulation(initial_value = .x, 
                                       N             = 240, 
                                       mean_return   = mean_port_return, 
                                       sd_return     = stddev_port_return)) %>%
    
    # Add column worth
    mutate(month = 1:nrow(.)) %>%
    select(month, everything()) %>%
    
    # Rerange column names
    set_names(c("month", names(starts))) %>%
    
    # Transform to long form
    pivot_longer(cols = -month, names_to = "sim", values_to = "growth")

Portfolio_sim
## # A tibble: 12,291 × 3
##    month sim   growth
##    <int> <chr>  <dbl>
##  1     1 sim1       1
##  2     1 sim2       1
##  3     1 sim3       1
##  4     1 sim4       1
##  5     1 sim5       1
##  6     1 sim6       1
##  7     1 sim7       1
##  8     1 sim8       1
##  9     1 sim9       1
## 10     1 sim10      1
## # ℹ 12,281 more rows
Portfolio_sim %>%
    
    group_by(sim) %>%
    summarize(growth = last(growth)) %>%
    ungroup() %>%
    pull(growth) %>%
    
    quantile(probs = c(0, 0.25, 0.5, 0.75, 1)) %>%
    round(2)
##       0%      25%      50%      75%     100% 
##    25.66   243.14   849.18  2076.95 20046.54

8 Visualizing simulations with ggplot

Line Plot of Simulations with Max, Median, and Min

sim_summary <- Portfolio_sim %>%
    
    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 20047.   849.  25.7
# Step 2
Portfolio_sim %>%
    
    # Filter by max, median, min
    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 1s over 20 months",
         subtitle = "Maximum, Median, and Minimum Simulation")
## $title
## [1] "Simulating growth of 1s over 20 months"
## 
## $subtitle
## [1] "Maximum, Median, and Minimum Simulation"
## 
## attr(,"class")
## [1] "labels"

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? ## Answer According to the Monte Carlo simulation results for a 100 investment over 20 years, the projected return is roughly 5,000, as illustrated by the green line. The red line represents the best-case scenario, which estimates the investment may reach 7500. The worst-case scenario, shown by the blue line, indicates little to no growth, with the investment remaining close to the original 100. However, this methodology is limited by its dependence on historical data and static inputs, which may not correctly reflect future market circumstances or capture severe volatility and unusual occurrences.