# Load necessary libraries
library(quantmod)
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(quadprog)
library(Matrix)
library(PortfolioAnalytics)
## Loading required package: foreach
## Loading required package: PerformanceAnalytics
##
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
##
## legend
## Registered S3 method overwritten by 'PortfolioAnalytics':
## method from
## print.constraint ROI
library(zoo)
library(xts)
library(portfolioBacktest)
# Define stock tickers
stks <- c("AAPL", "MSFT", "AMZN", "NVDA", "GOOGL", "GOOG", "META", "BRK-B",
"TSLA", "UNH", "XOM", "JNJ", "JPM", "V", "PG", "LLY", "AVGO",
"MA", "HD", "MRK")
# Download stock data
for (sym in stks) {
if (sym == "BRK-B") {
getSymbols("BRK-B", src = "yahoo", from = "2019-01-01", to = "2023-05-31")
assign("BRK.B", `BRK-B`) # Rename it in the workspace
} else {
getSymbols(sym, src = "yahoo", from = "2019-01-01", to = "2023-05-31")
}
}
# Extract adjusted prices and ensure proper column names
prices <- do.call(merge, lapply(stks, function(sym) {
if (sym == "BRK-B") sym <- "BRK.B" # Use the renamed variable
tmp <- Ad(get(sym))
colnames(tmp) <- sym # Ensure each column has a proper name
return(tmp)
}))
# Handle missing values
prices <- na.locf(prices)
# Validate structure
print(dim(prices)) # Check dimensions
## [1] 1110 20
print(colnames(prices)) # Ensure all stocks are correctly named
## [1] "AAPL" "MSFT" "AMZN" "NVDA" "GOOGL" "GOOG" "META" "BRK.B" "TSLA"
## [10] "UNH" "XOM" "JNJ" "JPM" "V" "PG" "LLY" "AVGO" "MA"
## [19] "HD" "MRK"
print(head(prices)) # Check for missing or incorrect values
## AAPL MSFT AMZN NVDA GOOGL GOOG META BRK.B
## 2019-01-02 37.66718 95.11981 76.9565 3.377779 52.48309 52.04599 135.0435 202.80
## 2019-01-03 33.91525 91.62054 75.0140 3.173704 51.02953 50.56351 131.1220 191.66
## 2019-01-04 35.36307 95.88174 78.7695 3.377035 53.64701 53.28313 137.3029 195.20
## 2019-01-07 35.28437 96.00403 81.4755 3.555819 53.54003 53.16767 137.4024 196.91
## 2019-01-08 35.95699 96.70012 82.8290 3.467295 54.01028 53.56031 141.8614 196.31
## 2019-01-09 36.56760 98.08289 82.9710 3.535485 53.82516 53.47970 143.5534 196.37
## TSLA UNH XOM JNJ JPM V PG
## 2019-01-02 20.67467 221.5286 51.71987 107.2963 82.91622 127.2749 77.99895
## 2019-01-03 20.02400 215.4874 50.92578 105.5913 81.73783 122.6883 77.45208
## 2019-01-04 21.17933 218.0076 52.80340 107.3635 84.75113 127.9739 79.03291
## 2019-01-07 22.33067 218.4261 53.07800 106.6748 84.81005 130.2816 78.71674
## 2019-01-08 22.35667 221.3466 53.46391 109.1525 84.65011 130.9902 79.00729
## 2019-01-09 22.56867 221.6650 53.74591 108.2874 84.50704 132.5318 77.71698
## LLY AVGO MA HD MRK
## 2019-01-02 105.1144 21.12924 183.2167 147.8488 59.73201
## 2019-01-03 101.8479 19.24977 174.9510 144.5901 58.50718
## 2019-01-04 104.9131 19.43897 183.2360 148.8865 60.26936
## 2019-01-07 105.4804 19.83486 184.6459 151.8193 59.60558
## 2019-01-08 106.4503 19.67567 185.9904 152.5481 60.04810
## 2019-01-09 107.2464 20.52665 189.3662 154.1260 59.58979
# Ensure no NA values before proceeding
if (any(is.na(prices))) {
stop("Error: 'prices' contains NA values. Fix before running backtest.")
}
# Save data
save(prices, file = "STKS_riskParity.RData")
# Load saved data
load("STKS_riskParity.RData")
# Risk parity function
risk_parity <- function(dataset) {
log_returns <- diff(log(dataset))[-1]
log_returns[is.na(log_returns)] <- 0
riskParityPortfolio(cov(log_returns))$w
}
# Max Sharpe ratio function
max_sharpe_ratio <- function(dataset) {
log_returns <- diff(log(dataset))[-1]
log_returns[is.na(log_returns)] <- 0
N <- ncol(dataset)
Sigma <- cov(log_returns)
mu <- colMeans(log_returns)
if (all(mu <= 1e-8)) return(rep(0, N))
Dmat <- 2 * Sigma
Amat <- cbind(mu, diag(N))
bvec <- c(1, rep(0, N))
dvec <- rep(0, N)
res <- try(solve.QP(Dmat = Dmat, dvec = dvec, Amat = Amat, bvec = bvec, meq = 1), silent = TRUE)
if (inherits(res, 'try-error')) {
D_mat <- nearPD(Dmat)
res <- solve.QP(Dmat = as.matrix(D_mat$mat), dvec = dvec, Amat = Amat, bvec = bvec, meq = 1)
}
w <- res$solution
w / sum(w)
}
# Wrap prices in a properly formatted list
stk_data <- list(adjusted = prices)
# Backtesting
bt <- portfolioBacktest(
list("Risk Parity" = risk_parity, "Tangency Portfolio" = max_sharpe_ratio),
list(stk_data), # Ensure dataset is a list of lists
lookback = 6,
optimize_every = 3,
rebalance_every = 1,
shortselling = TRUE,
benchmarks = '1/N',
execution = "next period"
)
## Backtesting 2 portfolios over 1 datasets (periodicity = daily data)...
## Backtesting benchmarks...
# Debug structure of backtest results
print(str(bt))
## List of 3
## $ Risk Parity :List of 1
## ..$ data1:List of 4
## .. ..$ performance : Named logi [1:13] NA NA NA NA NA NA ...
## .. .. ..- attr(*, "names")= chr [1:13] "Sharpe ratio" "max drawdown" "annual return" "annual volatility" ...
## .. .. ..- attr(*, "desired_direction")= num [1:13] 1 -1 1 -1 1 -1 1 1 -1 -1 ...
## .. ..$ cpu_time : logi NA
## .. ..$ error : logi TRUE
## .. ..$ error_message: chr "unused argument (w_current = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))"
## .. .. ..- attr(*, "error_stack")=List of 2
## .. .. .. ..$ at : chr [1:44] "(function (dataset) " "{" " log_returns <- diff(log(dataset))[-1]" " log_returns[is.na(log_returns)] <- 0" ...
## .. .. .. ..$ stack: chr ""
## $ Tangency Portfolio:List of 1
## ..$ data1:List of 4
## .. ..$ performance : Named logi [1:13] NA NA NA NA NA NA ...
## .. .. ..- attr(*, "names")= chr [1:13] "Sharpe ratio" "max drawdown" "annual return" "annual volatility" ...
## .. .. ..- attr(*, "desired_direction")= num [1:13] 1 -1 1 -1 1 -1 1 1 -1 -1 ...
## .. ..$ cpu_time : logi NA
## .. ..$ error : logi TRUE
## .. ..$ error_message: chr "unused argument (w_current = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))"
## .. .. ..- attr(*, "error_stack")=List of 2
## .. .. .. ..$ at : chr [1:61] "(function (dataset) " "{" " log_returns <- diff(log(dataset))[-1]" " log_returns[is.na(log_returns)] <- 0" ...
## .. .. .. ..$ stack: chr ""
## $ 1/N :List of 1
## ..$ data1:List of 10
## .. ..$ performance : Named num [1:13] 1.143 0.318 0.271 0.237 1.627 ...
## .. .. ..- attr(*, "names")= chr [1:13] "Sharpe ratio" "max drawdown" "annual return" "annual volatility" ...
## .. .. ..- attr(*, "desired_direction")= num [1:13] 1 -1 1 -1 1 -1 1 1 -1 -1 ...
## .. ..$ cpu_time : num 0.00101
## .. ..$ error : logi FALSE
## .. ..$ error_message: logi NA
## .. ..$ w_optimized :An xts object on 2019-01-09 / 2023-05-24 containing:
## Data: double [368, 20]
## Columns: AAPL, MSFT, AMZN, NVDA, GOOGL ... with 15 more columns
## Index: Date [368] (TZ: "UTC")
## .. ..$ w_rebalanced :An xts object on 2019-01-09 / 2023-05-25 containing:
## Data: double [1103, 20]
## Columns: AAPL, MSFT, AMZN, NVDA, GOOGL ... with 15 more columns
## Index: Date [1103] (TZ: "UTC")
## .. ..$ w_bop :An xts object on 2019-01-11 / 2023-05-30 containing:
## Data: double [1103, 20]
## Columns: AAPL, MSFT, AMZN, NVDA, GOOGL ... with 15 more columns
## Index: Date [1103] (TZ: "UTC")
## .. ..$ return :An xts object on 2019-01-11 / 2023-05-30 containing:
## Data: double [1103, 1]
## Columns: return
## Index: Date [1103] (TZ: "UTC")
## xts Attributes:
## .. .. .. $ ret_type : chr "discrete"
## .. .. .. $ coredata_content: chr "discreteReturn"
## .. ..$ wealth :An xts object on 2019-01-10 / 2023-05-30 containing:
## Data: double [1104, 1]
## Columns: NAV
## Index: Date [1104] (TZ: "UTC")
## .. ..$ X_lin :An xts object on 2019-01-02 / 2023-05-30 containing:
## Data: double [1110, 20]
## Columns: AAPL, MSFT, AMZN, NVDA, GOOGL ... with 15 more columns
## Index: Date [1110] (TZ: "UTC")
## xts Attributes:
## .. .. .. $ src : chr "yahoo"
## .. .. .. $ updated : POSIXct[1:1], format: "2025-04-01 07:37:59"
## .. .. .. $ ret_type : chr "discrete"
## .. .. .. $ coredata_content: chr "discreteReturn"
## - attr(*, "contain_benchmark")= logi TRUE
## - attr(*, "benchmark_index")= int 3
## NULL
# Validate backtest results before plotting
if (!is.null(bt$`Risk Parity`$data1$cum_return)) {
portfolioBacktest::backtestChartCumReturn(bt)
portfolioBacktest::backtestChartDrawdown(bt)
} else {
cat("Error: Backtest results are not properly formatted. Check input data.\n")
}
## Error: Backtest results are not properly formatted. Check input data.
# Weight Allocation Over Time
if (!is.null(bt$`Risk Parity`$data1$w_rebalanced)) {
backtestChartStackedBar(bt, portfolio = "Risk Parity", legend = TRUE)
backtestChartStackedBar(bt, portfolio = "Tangency Portfolio", legend = TRUE)
} else {
cat("Error: No rebalanced weights found in backtest.\n")
}
## Error: No rebalanced weights found in backtest.