#Accessing the keys
alpha_key <- Sys.getenv("ALPHAVANTAGE_KEY")
fred_key <- Sys.getenv("FRED_KEY")mp04
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
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.2library(dplyr)
Attaching package: 'dplyr'The following objects are masked from 'package:stats':
filter, lagThe 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. 26906629FRED 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 columnWarning: 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.567ggplot(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 forwardlibrary(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.1232Visualizing 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.01572207us_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.01970415Visualize 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.1269836Monte 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.703Overlay 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.0454516938Calculate 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 0Visualize 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:
Probability of Exhausting Savings:
- The first chart shows 0 occurrences of exhausted savings (all
FALSEforexhausted), suggesting that with the assumptions and bootstrap histories used, none of the retirement scenarios result in portfolio depletion.
- The first chart shows 0 occurrences of exhausted savings (all
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
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.
Starting Salary:
Higher salaries favor ORP for greater investment contributions.
Lower salaries may make TRS more attractive for guaranteed income.
Risk Tolerance:
TRS suits risk-averse individuals.
ORP suits those willing to accept variability for potential upside.
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
Bootstrap Sampling:
Historical data used in simulations may not reflect future economic conditions.
Bootstrap methods assume past trends are indicative of future outcomes.
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.
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.