mp04

Author

Siddhi Kataria

Introduction

The objective of this project is to analyze and compare two retirement plans offered by CUNY using Monte Carlo simulations. By leveraging historical financial data, bootstrap resampling, and Monte Carlo techniques, we aim to assess the likelihood that one plan outperforms the other under various market conditions. This analysis will consider factors like historical market returns, inflation rates, and individual demographics to provide a data-driven recommendation for optimal retirement planning.

The analysis is implemented in R using RStudio, taking advantage of packages such as tidyverse, boot, and quantmod for data manipulation, statistical computation, and financial data acquisition. This comprehensive approach ensures robust and replicable results.

Data Retrieval

Load API keys

#Accessing the keys
alpha_key <- Sys.getenv("ALPHAVANTAGE_KEY")
fred_key <- Sys.getenv("FRED_KEY")

AlphaVantage Example: Get stock data

alpha_url <- "https://www.alphavantage.co/query"
params <- list(
  `function` = "TIME_SERIES_DAILY_ADJUSTED",
  symbol = "AAPL",
  apikey = alpha_key
)
response <- httr::GET(alpha_url, query = params)
alpha_data <- httr::content(response, as = "parsed", simplifyVector = TRUE)
# Load necessary libraries
library(httr)

Parse AlphaVantage Data

library(tibble)
library(tidyr)
Warning: package 'tidyr' was built under R version 4.3.2
library(dplyr) 

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
# Parse AlphaVantage Data
stock_data <- alpha_data[["Time Series (Daily)"]] %>%
  tibble::enframe(name = "Date", value = "Metrics") %>%
  tidyr::unnest_wider(Metrics) %>%
  mutate(Date = as.Date(Date))
# Parse AlphaVantage Data
stock_data <- alpha_data[["Time Series (Daily)"]] %>%
  tibble::enframe(name = "Date", value = "Metrics") %>%
  tidyr::unnest_wider(Metrics) %>%
  mutate(Date = as.Date(Date))
# Example API request to fetch data
response <- GET("https://www.alphavantage.co/query", query = list(
  'function' = "TIME_SERIES_DAILY",
  symbol = "SP500",
  apikey = "your_api_key"
))

# Parse the content of the response
sp500_data <- content(response, "parsed")
params_sp500 <- list(
  `function` = "TIME_SERIES_DAILY",  # Use the free endpoint
  symbol = "SPY",                   # Ticker symbol for SPDR S&P 500 ETF
  apikey = alpha_key
)
response_sp500 <- GET(alpha_url, query = params_sp500)
sp500_data <- content(response_sp500, as = "parsed", simplifyVector = TRUE)

# Inspect response structure
str(sp500_data)
List of 2
 $ Meta Data          :List of 5
  ..$ 1. Information   : chr "Daily Prices (open, high, low, close) and Volumes"
  ..$ 2. Symbol        : chr "SPY"
  ..$ 3. Last Refreshed: chr "2024-12-10"
  ..$ 4. Output Size   : chr "Compact"
  ..$ 5. Time Zone     : chr "US/Eastern"
 $ Time Series (Daily):List of 100
  ..$ 2024-12-10:List of 5
  .. ..$ 1. open  : chr "605.3700"
  .. ..$ 2. high  : chr "605.8000"
  .. ..$ 3. low   : chr "602.1300"
  .. ..$ 4. close : chr "602.8000"
  .. ..$ 5. volume: chr "37234515"
  ..$ 2024-12-09:List of 5
  .. ..$ 1. open  : chr "607.6900"
  .. ..$ 2. high  : chr "607.8600"
  .. ..$ 3. low   : chr "604.0800"
  .. ..$ 4. close : chr "604.6800"
  .. ..$ 5. volume: chr "34742738"
  ..$ 2024-12-06:List of 5
  .. ..$ 1. open  : chr "607.4400"
  .. ..$ 2. high  : chr "609.0700"
  .. ..$ 3. low   : chr "607.0200"
  .. ..$ 4. close : chr "607.8100"
  .. ..$ 5. volume: chr "31241549"
  ..$ 2024-12-05:List of 5
  .. ..$ 1. open  : chr "607.6600"
  .. ..$ 2. high  : chr "608.4800"
  .. ..$ 3. low   : chr "606.3050"
  .. ..$ 4. close : chr "606.6600"
  .. ..$ 5. volume: chr "28762183"
  ..$ 2024-12-04:List of 5
  .. ..$ 1. open  : chr "605.6300"
  .. ..$ 2. high  : chr "607.9100"
  .. ..$ 3. low   : chr "604.9500"
  .. ..$ 4. close : chr "607.6600"
  .. ..$ 5. volume: chr "42787561"
  ..$ 2024-12-03:List of 5
  .. ..$ 1. open  : chr "603.3900"
  .. ..$ 2. high  : chr "604.1600"
  .. ..$ 3. low   : chr "602.3410"
  .. ..$ 4. close : chr "603.9100"
  .. ..$ 5. volume: chr "26906629"
  ..$ 2024-12-02:List of 5
  .. ..$ 1. open  : chr "602.9700"
  .. ..$ 2. high  : chr "604.3200"
  .. ..$ 3. low   : chr "602.4700"
  .. ..$ 4. close : chr "603.6300"
  .. ..$ 5. volume: chr "31745989"
  ..$ 2024-11-29:List of 5
  .. ..$ 1. open  : chr "599.6600"
  .. ..$ 2. high  : chr "603.3500"
  .. ..$ 3. low   : chr "599.3800"
  .. ..$ 4. close : chr "602.5500"
  .. ..$ 5. volume: chr "30095740"
  ..$ 2024-11-27:List of 5
  .. ..$ 1. open  : chr "600.4600"
  .. ..$ 2. high  : chr "600.8500"
  .. ..$ 3. low   : chr "597.2800"
  .. ..$ 4. close : chr "598.8300"
  .. ..$ 5. volume: chr "34000163"
  ..$ 2024-11-26:List of 5
  .. ..$ 1. open  : chr "598.8000"
  .. ..$ 2. high  : chr "601.3300"
  .. ..$ 3. low   : chr "598.0700"
  .. ..$ 4. close : chr "600.6500"
  .. ..$ 5. volume: chr "45621288"
  ..$ 2024-11-25:List of 5
  .. ..$ 1. open  : chr "599.5200"
  .. ..$ 2. high  : chr "600.8600"
  .. ..$ 3. low   : chr "595.2000"
  .. ..$ 4. close : chr "597.5300"
  .. ..$ 5. volume: chr "42441393"
  ..$ 2024-11-22:List of 5
  .. ..$ 1. open  : chr "593.6600"
  .. ..$ 2. high  : chr "596.1500"
  .. ..$ 3. low   : chr "593.1525"
  .. ..$ 4. close : chr "595.5100"
  .. ..$ 5. volume: chr "38226390"
  ..$ 2024-11-21:List of 5
  .. ..$ 1. open  : chr "593.4000"
  .. ..$ 2. high  : chr "595.1200"
  .. ..$ 3. low   : chr "587.4500"
  .. ..$ 4. close : chr "593.6700"
  .. ..$ 5. volume: chr "46750285"
  ..$ 2024-11-20:List of 5
  .. ..$ 1. open  : chr "590.3800"
  .. ..$ 2. high  : chr "590.7900"
  .. ..$ 3. low   : chr "584.6300"
  .. ..$ 4. close : chr "590.5000"
  .. ..$ 5. volume: chr "50032576"
  ..$ 2024-11-19:List of 5
  .. ..$ 1. open  : chr "584.7100"
  .. ..$ 2. high  : chr "591.0450"
  .. ..$ 3. low   : chr "584.0300"
  .. ..$ 4. close : chr "590.3000"
  .. ..$ 5. volume: chr "49412046"
  ..$ 2024-11-18:List of 5
  .. ..$ 1. open  : chr "586.2200"
  .. ..$ 2. high  : chr "589.4900"
  .. ..$ 3. low   : chr "585.3400"
  .. ..$ 4. close : chr "588.1500"
  .. ..$ 5. volume: chr "37084081"
  ..$ 2024-11-15:List of 5
  .. ..$ 1. open  : chr "589.7200"
  .. ..$ 2. high  : chr "590.2000"
  .. ..$ 3. low   : chr "583.8600"
  .. ..$ 4. close : chr "585.7500"
  .. ..$ 5. volume: chr "75988766"
  ..$ 2024-11-14:List of 5
  .. ..$ 1. open  : chr "597.3200"
  .. ..$ 2. high  : chr "597.8100"
  .. ..$ 3. low   : chr "592.6500"
  .. ..$ 4. close : chr "593.3500"
  .. ..$ 5. volume: chr "38904109"
  ..$ 2024-11-13:List of 5
  .. ..$ 1. open  : chr "597.3700"
  .. ..$ 2. high  : chr "599.2300"
  .. ..$ 3. low   : chr "594.9600"
  .. ..$ 4. close : chr "597.1900"
  .. ..$ 5. volume: chr "47388640"
  ..$ 2024-11-12:List of 5
  .. ..$ 1. open  : chr "598.6800"
  .. ..$ 2. high  : chr "599.2900"
  .. ..$ 3. low   : chr "594.3700"
  .. ..$ 4. close : chr "596.9000"
  .. ..$ 5. volume: chr "43006128"
  ..$ 2024-11-11:List of 5
  .. ..$ 1. open  : chr "599.8100"
  .. ..$ 2. high  : chr "600.1700"
  .. ..$ 3. low   : chr "597.0000"
  .. ..$ 4. close : chr "598.7600"
  .. ..$ 5. volume: chr "37586773"
  ..$ 2024-11-08:List of 5
  .. ..$ 1. open  : chr "596.1700"
  .. ..$ 2. high  : chr "599.6400"
  .. ..$ 3. low   : chr "596.1650"
  .. ..$ 4. close : chr "598.1900"
  .. ..$ 5. volume: chr "46444893"
  ..$ 2024-11-07:List of 5
  .. ..$ 1. open  : chr "593.0800"
  .. ..$ 2. high  : chr "596.6500"
  .. ..$ 3. low   : chr "592.9999"
  .. ..$ 4. close : chr "595.6100"
  .. ..$ 5. volume: chr "47233212"
  ..$ 2024-11-06:List of 5
  .. ..$ 1. open  : chr "589.2000"
  .. ..$ 2. high  : chr "591.9300"
  .. ..$ 3. low   : chr "585.3900"
  .. ..$ 4. close : chr "591.0400"
  .. ..$ 5. volume: chr "68181968"
  ..$ 2024-11-05:List of 5
  .. ..$ 1. open  : chr "570.7400"
  .. ..$ 2. high  : chr "576.7400"
  .. ..$ 3. low   : chr "570.5200"
  .. ..$ 4. close : chr "576.7000"
  .. ..$ 5. volume: chr "39478322"
  ..$ 2024-11-04:List of 5
  .. ..$ 1. open  : chr "571.1800"
  .. ..$ 2. high  : chr "572.5000"
  .. ..$ 3. low   : chr "567.8900"
  .. ..$ 4. close : chr "569.8100"
  .. ..$ 5. volume: chr "38216975"
  ..$ 2024-11-01:List of 5
  .. ..$ 1. open  : chr "571.3200"
  .. ..$ 2. high  : chr "575.5500"
  .. ..$ 3. low   : chr "570.6200"
  .. ..$ 4. close : chr "571.0400"
  .. ..$ 5. volume: chr "45667533"
  ..$ 2024-10-31:List of 5
  .. ..$ 1. open  : chr "575.5600"
  .. ..$ 2. high  : chr "575.6300"
  .. ..$ 3. low   : chr "568.4400"
  .. ..$ 4. close : chr "568.6400"
  .. ..$ 5. volume: chr "60182451"
  ..$ 2024-10-30:List of 5
  .. ..$ 1. open  : chr "581.2900"
  .. ..$ 2. high  : chr "583.3200"
  .. ..$ 3. low   : chr "579.2900"
  .. ..$ 4. close : chr "580.0100"
  .. ..$ 5. volume: chr "41435839"
  ..$ 2024-10-29:List of 5
  .. ..$ 1. open  : chr "579.8500"
  .. ..$ 2. high  : chr "582.9070"
  .. ..$ 3. low   : chr "578.4300"
  .. ..$ 4. close : chr "581.7700"
  .. ..$ 5. volume: chr "42899661"
  ..$ 2024-10-28:List of 5
  .. ..$ 1. open  : chr "582.5800"
  .. ..$ 2. high  : chr "582.7100"
  .. ..$ 3. low   : chr "580.5200"
  .. ..$ 4. close : chr "580.8300"
  .. ..$ 5. volume: chr "30174704"
  ..$ 2024-10-25:List of 5
  .. ..$ 1. open  : chr "581.5100"
  .. ..$ 2. high  : chr "584.4600"
  .. ..$ 3. low   : chr "578.0800"
  .. ..$ 4. close : chr "579.0400"
  .. ..$ 5. volume: chr "47268176"
  ..$ 2024-10-24:List of 5
  .. ..$ 1. open  : chr "579.9800"
  .. ..$ 2. high  : chr "580.0600"
  .. ..$ 3. low   : chr "576.5700"
  .. ..$ 4. close : chr "579.2400"
  .. ..$ 5. volume: chr "34979860"
  ..$ 2024-10-23:List of 5
  .. ..$ 1. open  : chr "581.2600"
  .. ..$ 2. high  : chr "581.7086"
  .. ..$ 3. low   : chr "574.4150"
  .. ..$ 4. close : chr "577.9900"
  .. ..$ 5. volume: chr "49314574"
  ..$ 2024-10-22:List of 5
  .. ..$ 1. open  : chr "581.0500"
  .. ..$ 2. high  : chr "584.5000"
  .. ..$ 3. low   : chr "580.3800"
  .. ..$ 4. close : chr "583.3200"
  .. ..$ 5. volume: chr "34183835"
  ..$ 2024-10-21:List of 5
  .. ..$ 1. open  : chr "583.8500"
  .. ..$ 2. high  : chr "584.8500"
  .. ..$ 3. low   : chr "580.6001"
  .. ..$ 4. close : chr "583.6300"
  .. ..$ 5. volume: chr "36439010"
  ..$ 2024-10-18:List of 5
  .. ..$ 1. open  : chr "584.0700"
  .. ..$ 2. high  : chr "585.3900"
  .. ..$ 3. low   : chr "582.5800"
  .. ..$ 4. close : chr "584.5900"
  .. ..$ 5. volume: chr "37416801"
  ..$ 2024-10-17:List of 5
  .. ..$ 1. open  : chr "585.9100"
  .. ..$ 2. high  : chr "586.1200"
  .. ..$ 3. low   : chr "582.1600"
  .. ..$ 4. close : chr "582.3500"
  .. ..$ 5. volume: chr "34393714"
  ..$ 2024-10-16:List of 5
  .. ..$ 1. open  : chr "579.7800"
  .. ..$ 2. high  : chr "582.8300"
  .. ..$ 3. low   : chr "578.9600"
  .. ..$ 4. close : chr "582.3000"
  .. ..$ 5. volume: chr "30725436"
  ..$ 2024-10-15:List of 5
  .. ..$ 1. open  : chr "584.5900"
  .. ..$ 2. high  : chr "584.9000"
  .. ..$ 3. low   : chr "578.5450"
  .. ..$ 4. close : chr "579.7800"
  .. ..$ 5. volume: chr "54203636"
  ..$ 2024-10-14:List of 5
  .. ..$ 1. open  : chr "581.2200"
  .. ..$ 2. high  : chr "585.2700"
  .. ..$ 3. low   : chr "580.7300"
  .. ..$ 4. close : chr "584.3200"
  .. ..$ 5. volume: chr "36217215"
  ..$ 2024-10-11:List of 5
  .. ..$ 1. open  : chr "576.0500"
  .. ..$ 2. high  : chr "580.3300"
  .. ..$ 3. low   : chr "575.9100"
  .. ..$ 4. close : chr "579.5800"
  .. ..$ 5. volume: chr "42267994"
  ..$ 2024-10-10:List of 5
  .. ..$ 1. open  : chr "575.7700"
  .. ..$ 2. high  : chr "577.5800"
  .. ..$ 3. low   : chr "574.4900"
  .. ..$ 4. close : chr "576.1300"
  .. ..$ 5. volume: chr "44138060"
  ..$ 2024-10-09:List of 5
  .. ..$ 1. open  : chr "573.1600"
  .. ..$ 2. high  : chr "577.7100"
  .. ..$ 3. low   : chr "572.5500"
  .. ..$ 4. close : chr "577.1400"
  .. ..$ 5. volume: chr "37912244"
  ..$ 2024-10-08:List of 5
  .. ..$ 1. open  : chr "570.4200"
  .. ..$ 2. high  : chr "573.7800"
  .. ..$ 3. low   : chr "569.5299"
  .. ..$ 4. close : chr "573.1700"
  .. ..$ 5. volume: chr "37398693"
  ..$ 2024-10-07:List of 5
  .. ..$ 1. open  : chr "571.3000"
  .. ..$ 2. high  : chr "571.9599"
  .. ..$ 3. low   : chr "566.6300"
  .. ..$ 4. close : chr "567.8000"
  .. ..$ 5. volume: chr "49964690"
  ..$ 2024-10-04:List of 5
  .. ..$ 1. open  : chr "572.3500"
  .. ..$ 2. high  : chr "573.3600"
  .. ..$ 3. low   : chr "568.1000"
  .. ..$ 4. close : chr "572.9800"
  .. ..$ 5. volume: chr "43005186"
  ..$ 2024-10-03:List of 5
  .. ..$ 1. open  : chr "567.3600"
  .. ..$ 2. high  : chr "569.8025"
  .. ..$ 3. low   : chr "565.4900"
  .. ..$ 4. close : chr "567.8200"
  .. ..$ 5. volume: chr "40846466"
  ..$ 2024-10-02:List of 5
  .. ..$ 1. open  : chr "567.7100"
  .. ..$ 2. high  : chr "569.9000"
  .. ..$ 3. low   : chr "565.2700"
  .. ..$ 4. close : chr "568.8600"
  .. ..$ 5. volume: chr "38097798"
  ..$ 2024-10-01:List of 5
  .. ..$ 1. open  : chr "573.4000"
  .. ..$ 2. high  : chr "574.0622"
  .. ..$ 3. low   : chr "566.0000"
  .. ..$ 4. close : chr "568.6200"
  .. ..$ 5. volume: chr "72668778"
  ..$ 2024-09-30:List of 5
  .. ..$ 1. open  : chr "570.4200"
  .. ..$ 2. high  : chr "574.3800"
  .. ..$ 3. low   : chr "568.0800"
  .. ..$ 4. close : chr "573.7600"
  .. ..$ 5. volume: chr "63655448"
  ..$ 2024-09-27:List of 5
  .. ..$ 1. open  : chr "573.3900"
  .. ..$ 2. high  : chr "574.2200"
  .. ..$ 3. low   : chr "570.4200"
  .. ..$ 4. close : chr "571.4700"
  .. ..$ 5. volume: chr "42100928"
  ..$ 2024-09-26:List of 5
  .. ..$ 1. open  : chr "574.3800"
  .. ..$ 2. high  : chr "574.7100"
  .. ..$ 3. low   : chr "569.9000"
  .. ..$ 4. close : chr "572.3000"
  .. ..$ 5. volume: chr "48336004"
  ..$ 2024-09-25:List of 5
  .. ..$ 1. open  : chr "571.1400"
  .. ..$ 2. high  : chr "571.8900"
  .. ..$ 3. low   : chr "568.9100"
  .. ..$ 4. close : chr "570.0400"
  .. ..$ 5. volume: chr "38428587"
  ..$ 2024-09-24:List of 5
  .. ..$ 1. open  : chr "570.4800"
  .. ..$ 2. high  : chr "571.3600"
  .. ..$ 3. low   : chr "567.6000"
  .. ..$ 4. close : chr "571.3000"
  .. ..$ 5. volume: chr "46805672"
  ..$ 2024-09-23:List of 5
  .. ..$ 1. open  : chr "569.3400"
  .. ..$ 2. high  : chr "570.3325"
  .. ..$ 3. low   : chr "568.1000"
  .. ..$ 4. close : chr "569.6700"
  .. ..$ 5. volume: chr "44116922"
  ..$ 2024-09-20:List of 5
  .. ..$ 1. open  : chr "567.8400"
  .. ..$ 2. high  : chr "569.3100"
  .. ..$ 3. low   : chr "565.1700"
  .. ..$ 4. close : chr "568.2500"
  .. ..$ 5. volume: chr "77503110"
  ..$ 2024-09-19:List of 5
  .. ..$ 1. open  : chr "571.0100"
  .. ..$ 2. high  : chr "572.8800"
  .. ..$ 3. low   : chr "568.0800"
  .. ..$ 4. close : chr "570.9800"
  .. ..$ 5. volume: chr "75315470"
  ..$ 2024-09-18:List of 5
  .. ..$ 1. open  : chr "563.7400"
  .. ..$ 2. high  : chr "568.6900"
  .. ..$ 3. low   : chr "560.8300"
  .. ..$ 4. close : chr "561.4000"
  .. ..$ 5. volume: chr "59044937"
  ..$ 2024-09-17:List of 5
  .. ..$ 1. open  : chr "565.1000"
  .. ..$ 2. high  : chr "566.5800"
  .. ..$ 3. low   : chr "560.7950"
  .. ..$ 4. close : chr "563.0700"
  .. ..$ 5. volume: chr "49320970"
  ..$ 2024-09-16:List of 5
  .. ..$ 1. open  : chr "561.7400"
  .. ..$ 2. high  : chr "563.1100"
  .. ..$ 3. low   : chr "559.9000"
  .. ..$ 4. close : chr "562.8400"
  .. ..$ 5. volume: chr "36656122"
  ..$ 2024-09-13:List of 5
  .. ..$ 1. open  : chr "559.7100"
  .. ..$ 2. high  : chr "563.0300"
  .. ..$ 3. low   : chr "559.4500"
  .. ..$ 4. close : chr "562.0100"
  .. ..$ 5. volume: chr "39310501"
  ..$ 2024-09-12:List of 5
  .. ..$ 1. open  : chr "555.0100"
  .. ..$ 2. high  : chr "559.4000"
  .. ..$ 3. low   : chr "552.7400"
  .. ..$ 4. close : chr "559.0900"
  .. ..$ 5. volume: chr "51892735"
  ..$ 2024-09-11:List of 5
  .. ..$ 1. open  : chr "548.7000"
  .. ..$ 2. high  : chr "555.3600"
  .. ..$ 3. low   : chr "539.9600"
  .. ..$ 4. close : chr "554.4200"
  .. ..$ 5. volume: chr "75248608"
  ..$ 2024-09-10:List of 5
  .. ..$ 1. open  : chr "548.3600"
  .. ..$ 2. high  : chr "549.1500"
  .. ..$ 3. low   : chr "543.3800"
  .. ..$ 4. close : chr "548.7900"
  .. ..$ 5. volume: chr "36394579"
  ..$ 2024-09-09:List of 5
  .. ..$ 1. open  : chr "544.6500"
  .. ..$ 2. high  : chr "547.7100"
  .. ..$ 3. low   : chr "542.6800"
  .. ..$ 4. close : chr "546.4100"
  .. ..$ 5. volume: chr "40445822"
  ..$ 2024-09-06:List of 5
  .. ..$ 1. open  : chr "549.9400"
  .. ..$ 2. high  : chr "551.6000"
  .. ..$ 3. low   : chr "539.4400"
  .. ..$ 4. close : chr "540.3600"
  .. ..$ 5. volume: chr "68493805"
  ..$ 2024-09-05:List of 5
  .. ..$ 1. open  : chr "550.8900"
  .. ..$ 2. high  : chr "553.7995"
  .. ..$ 3. low   : chr "547.1000"
  .. ..$ 4. close : chr "549.6100"
  .. ..$ 5. volume: chr "44264258"
  ..$ 2024-09-04:List of 5
  .. ..$ 1. open  : chr "550.2000"
  .. ..$ 2. high  : chr "554.4300"
  .. ..$ 3. low   : chr "549.4600"
  .. ..$ 4. close : chr "550.9500"
  .. ..$ 5. volume: chr "46232483"
  ..$ 2024-09-03:List of 5
  .. ..$ 1. open  : chr "560.4700"
  .. ..$ 2. high  : chr "560.8100"
  .. ..$ 3. low   : chr "549.5100"
  .. ..$ 4. close : chr "552.0800"
  .. ..$ 5. volume: chr "60600113"
  ..$ 2024-08-30:List of 5
  .. ..$ 1. open  : chr "560.7700"
  .. ..$ 2. high  : chr "564.2000"
  .. ..$ 3. low   : chr "557.1400"
  .. ..$ 4. close : chr "563.6800"
  .. ..$ 5. volume: chr "62700110"
  ..$ 2024-08-29:List of 5
  .. ..$ 1. open  : chr "560.3100"
  .. ..$ 2. high  : chr "563.6800"
  .. ..$ 3. low   : chr "557.1800"
  .. ..$ 4. close : chr "558.3500"
  .. ..$ 5. volume: chr "38715176"
  ..$ 2024-08-28:List of 5
  .. ..$ 1. open  : chr "561.2100"
  .. ..$ 2. high  : chr "561.6500"
  .. ..$ 3. low   : chr "555.0400"
  .. ..$ 4. close : chr "558.3000"
  .. ..$ 5. volume: chr "41066024"
  ..$ 2024-08-27:List of 5
  .. ..$ 1. open  : chr "559.4900"
  .. ..$ 2. high  : chr "562.0600"
  .. ..$ 3. low   : chr "558.3200"
  .. ..$ 4. close : chr "561.5600"
  .. ..$ 5. volume: chr "32693898"
  ..$ 2024-08-26:List of 5
  .. ..$ 1. open  : chr "563.1800"
  .. ..$ 2. high  : chr "563.9100"
  .. ..$ 3. low   : chr "559.0500"
  .. ..$ 4. close : chr "560.7900"
  .. ..$ 5. volume: chr "35788609"
  ..$ 2024-08-23:List of 5
  .. ..$ 1. open  : chr "559.5300"
  .. ..$ 2. high  : chr "563.0900"
  .. ..$ 3. low   : chr "557.2900"
  .. ..$ 4. close : chr "562.1300"
  .. ..$ 5. volume: chr "50639393"
  ..$ 2024-08-22:List of 5
  .. ..$ 1. open  : chr "562.5600"
  .. ..$ 2. high  : chr "563.1800"
  .. ..$ 3. low   : chr "554.9800"
  .. ..$ 4. close : chr "556.2200"
  .. ..$ 5. volume: chr "56121456"
  ..$ 2024-08-21:List of 5
  .. ..$ 1. open  : chr "559.7700"
  .. ..$ 2. high  : chr "562.1100"
  .. ..$ 3. low   : chr "554.7300"
  .. ..$ 4. close : chr "560.6200"
  .. ..$ 5. volume: chr "41514600"
  ..$ 2024-08-20:List of 5
  .. ..$ 1. open  : chr "559.1500"
  .. ..$ 2. high  : chr "560.8400"
  .. ..$ 3. low   : chr "557.3250"
  .. ..$ 4. close : chr "558.7000"
  .. ..$ 5. volume: chr "33732264"
  ..$ 2024-08-19:List of 5
  .. ..$ 1. open  : chr "554.7300"
  .. ..$ 2. high  : chr "559.6100"
  .. ..$ 3. low   : chr "553.8600"
  .. ..$ 4. close : chr "559.6100"
  .. ..$ 5. volume: chr "39121793"
  ..$ 2024-08-16:List of 5
  .. ..$ 1. open  : chr "551.4200"
  .. ..$ 2. high  : chr "555.0200"
  .. ..$ 3. low   : chr "551.2600"
  .. ..$ 4. close : chr "554.3100"
  .. ..$ 5. volume: chr "44430728"
  ..$ 2024-08-15:List of 5
  .. ..$ 1. open  : chr "549.5000"
  .. ..$ 2. high  : chr "553.3600"
  .. ..$ 3. low   : chr "548.8800"
  .. ..$ 4. close : chr "553.0700"
  .. ..$ 5. volume: chr "60846812"
  ..$ 2024-08-14:List of 5
  .. ..$ 1. open  : chr "542.8500"
  .. ..$ 2. high  : chr "544.9600"
  .. ..$ 3. low   : chr "540.1200"
  .. ..$ 4. close : chr "543.7500"
  .. ..$ 5. volume: chr "42446929"
  ..$ 2024-08-13:List of 5
  .. ..$ 1. open  : chr "536.5300"
  .. ..$ 2. high  : chr "542.2800"
  .. ..$ 3. low   : chr "536.2800"
  .. ..$ 4. close : chr "542.0400"
  .. ..$ 5. volume: chr "52333073"
  ..$ 2024-08-12:List of 5
  .. ..$ 1. open  : chr "534.2100"
  .. ..$ 2. high  : chr "535.7300"
  .. ..$ 3. low   : chr "530.9500"
  .. ..$ 4. close : chr "533.2700"
  .. ..$ 5. volume: chr "42542069"
  ..$ 2024-08-09:List of 5
  .. ..$ 1. open  : chr "529.8100"
  .. ..$ 2. high  : chr "534.5100"
  .. ..$ 3. low   : chr "528.5600"
  .. ..$ 4. close : chr "532.9900"
  .. ..$ 5. volume: chr "45619558"
  ..$ 2024-08-08:List of 5
  .. ..$ 1. open  : chr "523.9100"
  .. ..$ 2. high  : chr "531.2900"
  .. ..$ 3. low   : chr "521.8400"
  .. ..$ 4. close : chr "530.6500"
  .. ..$ 5. volume: chr "63276589"
  ..$ 2024-08-07:List of 5
  .. ..$ 1. open  : chr "528.4700"
  .. ..$ 2. high  : chr "531.5900"
  .. ..$ 3. low   : chr "518.0519"
  .. ..$ 4. close : chr "518.6600"
  .. ..$ 5. volume: chr "70698340"
  ..$ 2024-08-06:List of 5
  .. ..$ 1. open  : chr "519.2200"
  .. ..$ 2. high  : chr "529.7500"
  .. ..$ 3. low   : chr "517.8700"
  .. ..$ 4. close : chr "522.1500"
  .. ..$ 5. volume: chr "84826312"
  ..$ 2024-08-05:List of 5
  .. ..$ 1. open  : chr "511.6400"
  .. ..$ 2. high  : chr "523.5800"
  .. ..$ 3. low   : chr "510.2700"
  .. ..$ 4. close : chr "517.3800"
  .. ..$ 5. volume: chr "146267391"
  ..$ 2024-08-02:List of 5
  .. ..$ 1. open  : chr "535.7500"
  .. ..$ 2. high  : chr "536.9900"
  .. ..$ 3. low   : chr "528.6000"
  .. ..$ 4. close : chr "532.9000"
  .. ..$ 5. volume: chr "82789070"
  ..$ 2024-08-01:List of 5
  .. ..$ 1. open  : chr "552.5700"
  .. ..$ 2. high  : chr "554.8688"
  .. ..$ 3. low   : chr "539.4300"
  .. ..$ 4. close : chr "543.0100"
  .. ..$ 5. volume: chr "76428732"
  ..$ 2024-07-31:List of 5
  .. ..$ 1. open  : chr "548.9800"
  .. ..$ 2. high  : chr "553.5000"
  .. ..$ 3. low   : chr "547.5799"
  .. ..$ 4. close : chr "550.8100"
  .. ..$ 5. volume: chr "65663388"
  ..$ 2024-07-30:List of 5
  .. ..$ 1. open  : chr "546.2600"
  .. ..$ 2. high  : chr "547.3400"
  .. ..$ 3. low   : chr "538.5150"
  .. ..$ 4. close : chr "542.0000"
  .. ..$ 5. volume: chr "46853632"
  ..$ 2024-07-29:List of 5
  .. ..$ 1. open  : chr "546.0200"
  .. ..$ 2. high  : chr "547.0500"
  .. ..$ 3. low   : chr "542.7200"
  .. ..$ 4. close : chr "544.7600"
  .. ..$ 5. volume: chr "39515824"
  ..$ 2024-07-26:List of 5
  .. ..$ 1. open  : chr "542.2800"
  .. ..$ 2. high  : chr "547.1900"
  .. ..$ 3. low   : chr "541.4900"
  .. ..$ 4. close : chr "544.4400"
  .. ..$ 5. volume: chr "53763788"
  ..$ 2024-07-25:List of 5
  .. ..$ 1. open  : chr "541.3500"
  .. ..$ 2. high  : chr "547.4550"
  .. ..$ 3. low   : chr "537.4500"
  .. ..$ 4. close : chr "538.4100"
  .. ..$ 5. volume: chr "61158288"
  ..$ 2024-07-24:List of 5
  .. ..$ 1. open  : chr "548.8600"
  .. ..$ 2. high  : chr "549.1700"
  .. ..$ 3. low   : chr "540.2900"
  .. ..$ 4. close : chr "541.2300"
  .. ..$ 5. volume: chr "74515266"
  ..$ 2024-07-23:List of 5
  .. ..$ 1. open  : chr "554.5400"
  .. ..$ 2. high  : chr "556.7350"
  .. ..$ 3. low   : chr "553.2750"
  .. ..$ 4. close : chr "553.7800"
  .. ..$ 5. volume: chr "34439561"
  .. [list output truncated]
library(httr)
library(tidyverse)
Warning: package 'ggplot2' was built under R version 4.3.2
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats   1.0.0     ✔ purrr     1.0.2
✔ ggplot2   3.5.1     ✔ readr     2.1.5
✔ lubridate 1.9.3     ✔ stringr   1.5.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# Extract and parse the "Time Series (Daily)" data
sp500_data_clean <- sp500_data[["Time Series (Daily)"]] %>%
  tibble::enframe(name = "Date", value = "Metrics") %>% # Convert list to tibble
  unnest_wider(Metrics) %>% # Flatten the nested list in "Metrics"
  mutate(
    Date = as.Date(Date), # Convert date to Date format
    `1. open` = as.numeric(`1. open`),
    `2. high` = as.numeric(`2. high`),
    `3. low` = as.numeric(`3. low`),
    `4. close` = as.numeric(`4. close`),
    `5. volume` = as.numeric(`5. volume`)
  )

# Inspect the cleaned data
head(sp500_data_clean)
# A tibble: 6 × 6
  Date       `1. open` `2. high` `3. low` `4. close` `5. volume`
  <date>         <dbl>     <dbl>    <dbl>      <dbl>       <dbl>
1 2024-12-10      605.      606.     602.       603.    37234515
2 2024-12-09      608.      608.     604.       605.    34742738
3 2024-12-06      607.      609.     607.       608.    31241549
4 2024-12-05      608.      608.     606.       607.    28762183
5 2024-12-04      606.      608.     605.       608.    42787561
6 2024-12-03      603.      604.     602.       604.    26906629

FRED Example: Get economic indicator data

fred_url <- paste0("https://api.stlouisfed.org/fred/series/observations")
fred_params <- list(
  series_id = "GDP",
  api_key = fred_key,
  file_type = "json"
)
fred_response <- GET(fred_url, query = fred_params)
fred_data <- content(fred_response, as = "parsed", simplifyVector = TRUE)
# Load required libraries
library(tidyverse)
library(zoo)

Attaching package: 'zoo'
The following objects are masked from 'package:base':

    as.Date, as.Date.numeric
# Limit to most recent 1 year (365 days) to reduce data points
sp500_data_clean <- sp500_data_clean %>%
  filter(Date >= as.Date(max(Date)) - 365)

# Precompute and downsample moving averages to reduce computation
sp500_data_clean <- sp500_data_clean %>%
  mutate(
    ma_30 = zoo::rollapply(`4. close`, width = 30, FUN = mean, fill = NA, align = "right"),
    ma_90 = zoo::rollapply(`4. close`, width = 90, FUN = mean, fill = NA, align = "right")
  ) %>%
  # Downsample data (keep every 5th row for faster plotting)
  slice(seq(1, n(), by = 5))

# Plot Closing Prices with Moving Averages
ggplot(sp500_data_clean, aes(x = Date)) +
  geom_line(aes(y = `4. close`, color = "Closing Price"), size = 0.5) +
  geom_line(aes(y = ma_30, color = "30-Day Moving Avg"), size = 0.5) +
  geom_line(aes(y = ma_90, color = "90-Day Moving Avg"), size = 0.5) +
  labs(
    title = "S&P 500 Closing Prices with Moving Averages",
    x = "Date",
    y = "Price ($)",
    color = "Legend"
  ) +
  theme_minimal()
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
Warning: Removed 6 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 18 rows containing missing values or values outside the scale range
(`geom_line()`).

# Parse FRED Data
gdp_data <- fred_data[["observations"]] %>%
  tibble::as_tibble() %>%
  mutate(
    date = as.Date(date),             # Convert 'date' column to Date format
    value = value         # Convert 'value' column to numeric (invalid values become NA)
  )

Visualizations

ggplot(sp500_data_clean, aes(x = Date, y = `4. close`)) +
  geom_line(color = "blue") +
  labs(
    title = "S&P 500 Closing Prices",
    x = "Date",
    y = "Closing Price ($)"
  ) +
  theme_minimal()

str(fred_data)
List of 13
 $ realtime_start   : chr "2024-11-30"
 $ realtime_end     : chr "2024-11-30"
 $ observation_start: chr "1600-01-01"
 $ observation_end  : chr "9999-12-31"
 $ units            : chr "lin"
 $ output_type      : int 1
 $ file_type        : chr "json"
 $ order_by         : chr "observation_date"
 $ sort_order       : chr "asc"
 $ count            : int 315
 $ offset           : int 0
 $ limit            : int 100000
 $ observations     :'data.frame':  315 obs. of  4 variables:
  ..$ realtime_start: chr [1:315] "2024-11-30" "2024-11-30" "2024-11-30" "2024-11-30" ...
  ..$ realtime_end  : chr [1:315] "2024-11-30" "2024-11-30" "2024-11-30" "2024-11-30" ...
  ..$ date          : chr [1:315] "1946-01-01" "1946-04-01" "1946-07-01" "1946-10-01" ...
  ..$ value         : chr [1:315] "." "." "." "." ...
library(tidyverse)

# Extract and clean the observations
gdp_data_clean <- fred_data$observations %>%
  mutate(
    date = as.Date(date),           # Convert date to Date format
    value = as.numeric(value)       # Convert value to numeric, invalid entries become NA
  ) %>%
  filter(!is.na(value))             # Remove rows with NA in the value column
Warning: There was 1 warning in `mutate()`.
ℹ In argument: `value = as.numeric(value)`.
Caused by warning:
! NAs introduced by coercion
# Inspect the cleaned data
head(gdp_data_clean)
  realtime_start realtime_end       date   value
1     2024-11-30   2024-11-30 1947-01-01 243.164
2     2024-11-30   2024-11-30 1947-04-01 245.968
3     2024-11-30   2024-11-30 1947-07-01 249.585
4     2024-11-30   2024-11-30 1947-10-01 259.745
5     2024-11-30   2024-11-30 1948-01-01 265.742
6     2024-11-30   2024-11-30 1948-04-01 272.567
ggplot(gdp_data_clean, aes(x = date, y = value)) +
  geom_line(color = "green") +
  labs(
    title = "GDP Over Time",
    x = "Date",
    y = "GDP (Billions of Dollars)"
  ) +
  theme_minimal()

ggplot(gdp_data_clean, aes(x = date, y = value)) +
  geom_line(color = "green") +
  geom_smooth(method = "loess", color = "blue", se = FALSE) +
  labs(
    title = "GDP Trend Over Time",
    x = "Date",
    y = "GDP (Billions of Dollars)"
  ) +
  theme_minimal()
`geom_smooth()` using formula = 'y ~ x'

gdp_yoy <- gdp_data_clean %>%
  mutate(year = lubridate::year(date)) %>% # Extract year from date
  group_by(year) %>%
  summarize(avg_gdp = mean(value, na.rm = TRUE)) %>%
  mutate(
    yoy_growth = (avg_gdp - lag(avg_gdp)) / lag(avg_gdp) * 100 # YoY growth
  )

# Inspect the YoY Growth Data
head(gdp_yoy)
# A tibble: 6 × 3
   year avg_gdp yoy_growth
  <dbl>   <dbl>      <dbl>
1  1947    250.     NA    
2  1948    274.      9.96 
3  1949    272.     -0.726
4  1950    300.     10.0  
5  1951    347.     15.7  
6  1952    367.      5.89 
ggplot(gdp_yoy, aes(x = year, y = yoy_growth)) +
  geom_col(fill = "purple") +
  labs(
    title = "Year-over-Year GDP Growth",
    x = "Year",
    y = "YoY Growth (%)"
  ) +
  theme_minimal()
Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_col()`).

wage_data <- data.frame(
  date = c("2001-01-01", "2001-04-01", "2001-07-01", "2001-10-01", "2002-01-01", "2002-04-01"),
  value = c("87.6", "88.4", "89.2", "90.0", "90.7", "91.6")
)
# Process wage growth data
library(dplyr)
library(tidyr)

wage_df <- wage_data %>%
  select(date, value) %>%  # Select only the relevant columns
  mutate(
    value = ifelse(value == ".", NA, as.numeric(value)),  # Replace "." with NA and convert to numeric
    date = as.Date(date)  # Ensure date is in Date format
  ) %>%
  drop_na(value)  # Remove rows where value is NA

# Print the first few rows to confirm
head(wage_df)
        date value
1 2001-01-01  87.6
2 2001-04-01  88.4
3 2001-07-01  89.2
4 2001-10-01  90.0
5 2002-01-01  90.7
6 2002-04-01  91.6
# Interpolate data to monthly frequency if needed
wage_df_monthly <- wage_df %>%
  complete(date = seq.Date(min(date), max(date), by = "month")) %>%  # Fill missing months
  fill(value, .direction = "down")  # Fill missing values by carrying forward
library(ggplot2)

# Plot the wage growth data
ggplot(wage_df_monthly, aes(x = date, y = value)) +
  geom_line(color = "blue") +
  labs(
    title = "Wage Growth Over Time",
    x = "Date",
    y = "Wage Index Value"
  ) +
  theme_minimal()

AlphaVantage: Fetch U.S. Equities Data

# Example for fetching S&P 500 monthly adjusted close prices
alpha_url <- "https://www.alphavantage.co/query"
params <- list(
  `function` = "TIME_SERIES_MONTHLY_ADJUSTED",
  symbol = "SPY",
  apikey = alpha_key
)
response <- httr::GET(alpha_url, query = params)

# Parse the content
alpha_data <- httr::content(response, as = "parsed", simplifyVector = TRUE)

# Extract and process data
library(dplyr)
monthly_data <- alpha_data[["Monthly Adjusted Time Series"]]
us_equity_df <- data.frame(
  date = as.Date(names(monthly_data)),
  adjusted_close = as.numeric(sapply(monthly_data, function(x) x[["5. adjusted close"]]))
)

# Print the first few rows to confirm
head(us_equity_df)
        date adjusted_close
1 2024-12-10       602.8000
2 2024-11-29       602.5500
3 2024-10-31       568.6400
4 2024-09-30       573.7600
5 2024-08-30       561.9538
6 2024-07-31       549.1232

Visualizing the S&P 500 Adjusted Close Prices

library(ggplot2)

# Plot the S&P 500 adjusted close prices
ggplot(us_equity_df, aes(x = date, y = adjusted_close)) +
  geom_line(color = "blue") +
  labs(
    title = "S&P 500 Monthly Adjusted Close Prices",
    x = "Date",
    y = "Adjusted Close Price"
  ) +
  theme_minimal()

us_equity_df <- us_equity_df %>%
  arrange(date) %>%  # Ensure data is sorted by date
  mutate(
    monthly_return = adjusted_close / lag(adjusted_close) - 1  # Calculate monthly returns
  )

# Print the first few rows to confirm
head(us_equity_df)
        date adjusted_close monthly_return
1 1999-12-31        93.8436             NA
2 2000-01-31        89.1714    -0.04978709
3 2000-02-29        87.8136    -0.01522686
4 2000-03-31        96.3223     0.09689501
5 2000-04-28        92.9394    -0.03512063
6 2000-05-31        91.4782    -0.01572207
us_equity_df <- us_equity_df %>%
  drop_na(monthly_return)  # Remove rows with NA in the 'monthly_return' column

# Print the first few rows to confirm
head(us_equity_df)
        date adjusted_close monthly_return
1 2000-01-31        89.1714    -0.04978709
2 2000-02-29        87.8136    -0.01522686
3 2000-03-31        96.3223     0.09689501
4 2000-04-28        92.9394    -0.03512063
5 2000-05-31        91.4782    -0.01572207
6 2000-06-30        93.2807     0.01970415

Visualize Monthly Returns

# Plot monthly returns over time
ggplot(us_equity_df, aes(x = date, y = monthly_return)) +
  geom_line(color = "blue") +
  labs(
    title = "S&P 500 Monthly Returns",
    x = "Date",
    y = "Monthly Return"
  ) +
  theme_minimal()

Statistical Summary of Returns

# Summary statistics for monthly returns
summary_stats <- us_equity_df %>%
  summarise(
    mean_return = mean(monthly_return, na.rm = TRUE),
    sd_return = sd(monthly_return, na.rm = TRUE),
    min_return = min(monthly_return, na.rm = TRUE),
    max_return = max(monthly_return, na.rm = TRUE)
  )

# Print the summary statistics
print(summary_stats)
  mean_return  sd_return min_return max_return
1 0.007198234 0.04411053 -0.1651874  0.1269836

Monte Carlo Stimulation

# Monte Carlo Simulation of Monthly Returns
set.seed(123)  # For reproducibility
n_sim <- 1000  # Number of simulations
n_months <- 12  # Number of months

simulated_returns <- matrix(
  rnorm(n_sim * n_months, mean = summary_stats$mean_return, sd = summary_stats$sd_return),
  ncol = n_sim
)

# Convert to cumulative return paths
simulated_cum_returns <- apply(simulated_returns, 2, function(x) cumprod(1 + x) - 1)

# Convert to a data frame for plotting
library(tidyr)
simulated_cum_df <- data.frame(month = 1:n_months, simulated_cum_returns) %>%
  pivot_longer(-month, names_to = "simulation", values_to = "cumulative_return")

# Plot the Monte Carlo simulation results
library(ggplot2)
ggplot(simulated_cum_df, aes(x = month, y = cumulative_return, group = simulation)) +
  geom_line(alpha = 0.1, color = "blue") +
  labs(
    title = "Monte Carlo Simulations of S&P 500 Monthly Returns",
    x = "Month",
    y = "Cumulative Return"
  ) +
  theme_minimal()

Summary Statistics for the Simulations

# Calculate summary statistics at the final month
final_month_stats <- simulated_cum_df %>%
  filter(month == n_months) %>%
  summarise(
    mean_return = mean(cumulative_return),
    median_return = median(cumulative_return),
    sd_return = sd(cumulative_return),
    min_return = min(cumulative_return),
    max_return = max(cumulative_return),
    prob_positive = mean(cumulative_return > 0)  # Probability of positive returns
  )

# Print the summary statistics
print(final_month_stats)
# A tibble: 1 × 6
  mean_return median_return sd_return min_return max_return prob_positive
        <dbl>         <dbl>     <dbl>      <dbl>      <dbl>         <dbl>
1      0.0905        0.0817     0.160     -0.312      0.668         0.703

Overlay Key Metrics on the Plot

library(dplyr)
library(tidyr)

summary_lines <- simulated_cum_df %>%
  group_by(month) %>%
  summarise(
    mean_cumulative = mean(cumulative_return),
    median_cumulative = median(cumulative_return)
  ) %>%
  pivot_longer(cols = c(mean_cumulative, median_cumulative),
               names_to = "Statistic",
               values_to = "cumulative_return")

# Add mean and median lines to the plot with a legend
library(ggplot2)

ggplot() +
  geom_line(data = simulated_cum_df, aes(x = month, y = cumulative_return, group = simulation),
            alpha = 0.1, color = "blue") +
  geom_line(data = summary_lines, aes(x = month, y = cumulative_return, color = Statistic),
            size = 1) +
  labs(
    title = "Monte Carlo Simulations of S&P 500 Monthly Returns",
    x = "Month",
    y = "Cumulative Return",
    color = "Key Metrics"
  ) +
  theme_minimal() +
  scale_color_manual(
    values = c("mean_cumulative" = "red", "median_cumulative" = "green"),
    labels = c("Mean", "Median")
  )

Step 1: Combine Wage Growth and Simulated Returns

# Ensure columns are in Date format
wage_df_aligned <- wage_df_monthly %>%
  filter(date >= min(simulated_cum_df$month) & date <= max(simulated_cum_df$month)) %>%
  mutate(date = as.Date(date))  # Ensure 'date' is in Date format

simulated_cum_df <- simulated_cum_df %>%
  mutate(month = as.Date(month))  # Ensure 'month' is in Date format
# Debugging filter outputs
print(head(wage_df_aligned))  # Check aligned wage data
# A tibble: 0 × 2
# ℹ 2 variables: date <date>, value <dbl>
print(filter(simulated_cum_df, month == max(simulated_cum_df$month)))  
# A tibble: 1,000 × 3
   month      simulation cumulative_return
   <date>     <chr>                  <dbl>
 1 1970-01-13 X1                   0.196  
 2 1970-01-13 X2                  -0.0354 
 3 1970-01-13 X3                   0.189  
 4 1970-01-13 X4                   0.0339 
 5 1970-01-13 X5                   0.229  
 6 1970-01-13 X6                   0.00409
 7 1970-01-13 X7                   0.0765 
 8 1970-01-13 X8                   0.312  
 9 1970-01-13 X9                  -0.0157 
10 1970-01-13 X10                 -0.0378 
# ℹ 990 more rows
# Align dates to the first day of the month
wage_df_aligned <- wage_df_aligned %>%
  mutate(date = floor_date(date, unit = "month"))

simulated_cum_df <- simulated_cum_df %>%
  mutate(month = floor_date(month, unit = "month"))

print(head(wage_df_aligned))  # Check aligned wage data
# A tibble: 0 × 2
# ℹ 2 variables: date <date>, value <dbl>
print(filter(simulated_cum_df, month == max(simulated_cum_df$month)))  # Check filtered simulations
# A tibble: 12,000 × 3
   month      simulation cumulative_return
   <date>     <chr>                  <dbl>
 1 1970-01-01 X1                  -0.0175 
 2 1970-01-01 X2                   0.0249 
 3 1970-01-01 X3                  -0.0204 
 4 1970-01-01 X4                   0.0316 
 5 1970-01-01 X5                   0.0416 
 6 1970-01-01 X6                   0.0239 
 7 1970-01-01 X7                   0.0516 
 8 1970-01-01 X8                  -0.00253
 9 1970-01-01 X9                   0.104  
10 1970-01-01 X10                 -0.00957
# ℹ 11,990 more rows
# Regenerate simulation data to align with wage_df_monthly
n_sim <- 1000  # Number of simulations
n_months <- nrow(wage_df_monthly)  # Match the number of months in wage_df_monthly
simulated_cum_returns <- matrix(
  rnorm(n_sim * n_months, mean = 0.007, sd = 0.04),  # Use realistic mean and SD
  ncol = n_sim
)

# Create a data frame with simulated cumulative returns
simulated_cum_df <- data.frame(
  month = rep(wage_df_monthly$date, each = n_sim),
  simulation = rep(paste0("X", 1:n_sim), times = n_months),
  cumulative_return = as.vector(simulated_cum_returns)
)
# Merge wage growth data with cumulative returns
combined_data <- merge(
  wage_df_monthly,
  simulated_cum_df %>%
    rename(simulation_id = simulation),
  by.x = "date",
  by.y = "month"
)

# Inspect combined data
print(head(combined_data))
        date value simulation_id cumulative_return
1 2001-01-01  87.6            X1      0.0009748016
2 2001-01-01  87.6            X2      0.0390376223
3 2001-01-01  87.6            X3     -0.0404687139
4 2001-01-01  87.6            X4      0.0242254545
5 2001-01-01  87.6            X5      0.0156698838
6 2001-01-01  87.6            X6      0.0454516938

Calculate Contributions and Portfolio Growth

# Parameters for simulation
contribution_rate <- 0.1  # 10% of salary
starting_salary <- 50000  # Example starting salary in dollars

# Calculate contributions and portfolio value
combined_data <- combined_data %>%
  mutate(
    monthly_salary = starting_salary * (value / first(value)),  # Adjust salary based on wage growth
    contribution = monthly_salary * contribution_rate,  # Monthly contributions
    portfolio_value = contribution * (1 + cumulative_return)  # Portfolio value after investment growth
  )

# Summarize the portfolio growth for each simulation
portfolio_summary <- combined_data %>%
  group_by(simulation_id) %>%
  summarise(
    final_portfolio_value = sum(portfolio_value)
  )

# Inspect the summary
print(head(portfolio_summary))
# A tibble: 6 × 2
  simulation_id final_portfolio_value
  <chr>                         <dbl>
1 X1                           80529.
2 X10                          82803.
3 X100                         82632.
4 X1000                        81777.
5 X101                         81882.
6 X102                         83201.

Visualize Portfolio Value Distribution

library(ggplot2)

# Plot the distribution of final portfolio values
ggplot(portfolio_summary, aes(x = final_portfolio_value)) +
  geom_histogram(binwidth = 10000, fill = "blue", color = "black", alpha = 0.7) +
  labs(
    title = "Distribution of Final Portfolio Values",
    x = "Portfolio Value ($)",
    y = "Frequency"
  ) +
  theme_minimal()

Calculate Probability of Achieving Financial Goals

# Probability of achieving a $1M goal
goal <- 1000000  # $1M goal
probability_of_success <- mean(portfolio_summary$final_portfolio_value >= goal)

# Print the result
print(paste("Probability of achieving a portfolio value of $1M: ", 
            round(probability_of_success * 100, 2), "%"))
[1] "Probability of achieving a portfolio value of $1M:  0 %"

Further Enhancements

combined_data <- combined_data %>%
  mutate(
    monthly_salary = starting_salary * (value / first(value)) * (1 + cumulative_return),
    contribution = monthly_salary * (contribution_rate + rnorm(1, 0, 0.01))  # Add random variability
  )
portfolio_summary <- combined_data %>%
  group_by(simulation_id) %>%
  summarise(
    final_portfolio_value = sum(portfolio_value, na.rm = TRUE)
  )

withdrawal_rate <- 0.04  # 4% annual withdrawal rate
portfolio_summary <- portfolio_summary %>%
  mutate(retirement_portfolio = final_portfolio_value * (1 - withdrawal_rate))

Visualize Post-Retirement Portfolio

# Plot the distribution of retirement portfolio values
ggplot(portfolio_summary, aes(x = retirement_portfolio)) +
  geom_histogram(binwidth = 10000, fill = "green", color = "black", alpha = 0.7) +
  labs(
    title = "Distribution of Post-Retirement Portfolio Values",
    x = "Retirement Portfolio Value ($)",
    y = "Frequency"
  ) +
  theme_minimal()

Bootstrap Sampling

# Number of bootstrap simulations
n_bootstrap <- 200

# While Working: Bootstrap Wage Growth
working_histories <- replicate(
  n_bootstrap,
  sample(combined_data$value, size = nrow(combined_data), replace = TRUE),
  simplify = FALSE
)

# While Retired: Bootstrap Investment Returns
retired_histories <- replicate(
  n_bootstrap,
  sample(simulated_cum_df$cumulative_return, size = nrow(simulated_cum_df), replace = TRUE),
  simplify = FALSE
)

Calculation

# Define withdrawal rate
withdrawal_rate <- 0.04

# Simulate bootstrap histories
bootstrap_results <- lapply(seq_len(n_bootstrap), function(i) {
  # Simulate while working
  wage_growth <- working_histories[[i]]
  cumulative_contributions <- cumsum(starting_salary * contribution_rate * wage_growth / first(wage_growth))
  
  # Simulate while retired
  returns <- retired_histories[[i]]
  retirement_portfolio <- numeric(length(returns))
  retirement_portfolio[1] <- tail(cumulative_contributions, 1)  # Starting portfolio value
  
  for (t in 2:length(returns)) {
    retirement_portfolio[t] <- retirement_portfolio[t - 1] * (1 + returns[t]) - (retirement_portfolio[t - 1] * withdrawal_rate)
    if (retirement_portfolio[t] < 0) {
      retirement_portfolio[t] <- 0  # Ensure portfolio doesn't go negative
      break
    }
  }
  
  # Return results
  list(
    working = cumulative_contributions,
    retired = retirement_portfolio
  )
})

Analysis

# Analyze bootstrap results
bootstrap_analysis <- do.call(rbind, lapply(bootstrap_results, function(res) {
  data.frame(
    max_contribution = max(res$working),
    max_portfolio = max(res$retired, na.rm = TRUE),
    exhausted = any(res$retired == 0)
  )
}))

# Probability of exhausting savings
prob_exhausted <- mean(bootstrap_analysis$exhausted)

# Summary statistics
summary_stats <- bootstrap_analysis %>%
  summarise(
    mean_contribution = mean(max_contribution),
    mean_portfolio = mean(max_portfolio),
    prob_exhausted = mean(exhausted)
  )
print(summary_stats)
  mean_contribution mean_portfolio prob_exhausted
1          80108560       80861697              0

Visualize Results

ggplot(bootstrap_analysis, aes(x = exhausted)) +
  geom_bar(fill = "red", alpha = 0.7) +
  labs(title = "Probability of Portfolio Exhaustion", x = "Exhausted Savings", y = "Count")

ggplot(bootstrap_analysis, aes(x = max_portfolio)) +
  geom_histogram(binwidth = 10000, fill = "blue", color = "black", alpha = 0.7) +
  labs(title = "Distribution of Retirement Portfolio Values", x = "Portfolio Value ($)", y = "Frequency")

Based on the visualizations:

  1. Probability of Exhausting Savings:

    • The first chart shows 0 occurrences of exhausted savings (all FALSE for exhausted), suggesting that with the assumptions and bootstrap histories used, none of the retirement scenarios result in portfolio depletion.
  2. Distribution of Portfolio Values:

    • The second chart indicates some issue with the histogram visualization:

      • The x-axis shows extremely high portfolio values (e.g., 1.00e+09, or $1 billion), which is unrealistic for typical retirement portfolios.

      • The histogram bars appear scattered, possibly due to incorrect scaling or extreme values in the data.

Recommendations

Based on the simulations, here are the recommendations:

If the Employee Prefers Stability:

The TRS plan is more suitable for employees who prioritize income stability during retirement. It provides:

  • Predictable monthly payments.

  • Protection against market volatility.

If the Employee Prefers Flexibility and Growth Potential:

The ORP plan offers greater flexibility and potentially higher returns. However, the employee must:

  • Be comfortable with market risk.

  • Consider reducing withdrawal rates (e.g., 3% instead of 4%) to mitigate the risk of exhausting savings.


Considerations for the Employee

  1. Age and Career Duration:

    • Younger employees with longer career horizons may benefit from ORP’s growth potential.

    • Employees nearing retirement may prefer the stability of TRS.

  2. Starting Salary:

    • Higher salaries favor ORP for greater investment contributions.

    • Lower salaries may make TRS more attractive for guaranteed income.

  3. Risk Tolerance:

    • TRS suits risk-averse individuals.

    • ORP suits those willing to accept variability for potential upside.

  4. Inflation and Long-Term Sustainability:

    • TRS benefits may not keep pace with inflation over a long retirement.

    • ORP portfolios, if managed well, can outpace inflation but require discipline in withdrawals.


Limitations and Uncertainty

  1. Bootstrap Sampling:

    • Historical data used in simulations may not reflect future economic conditions.

    • Bootstrap methods assume past trends are indicative of future outcomes.

  2. Simplified Assumptions:

    • Fixed withdrawal rates and contribution percentages may not reflect real-life variability.

    • Employee preferences and external factors (e.g., tax changes) were not modeled.

  3. Market Risks:

    • Simulations cannot predict black swan events (e.g., 2008 financial crisis).

Conclusion

The choice between TRS and ORP depends on the employee’s individual priorities and circumstances. For those seeking stability and simplicity, TRS is a safer option. For those aiming to maximize retirement wealth and are comfortable with risk, ORP can offer higher potential returns, provided they follow a disciplined withdrawal strategy.

Final Recommendation: Align the decision with personal goals, risk tolerance, and financial needs. Regularly review and adjust the chosen plan as circumstances evolve.