I’m trying to find the portfolio with the maximum return at a given risk level using the fPortfolio package in R. However, there appears to be a bug in the package that the optimal portfolio (based on max returns) does not change when changing the target risk level, suggesting that it is not taking into account the target risk level specified by the user. This issue has been identified by others (see links below), and I provide a minimal reproducible example below. As far as I can tell, this issue has not been addressed in the package documentation (archived here) or in the authors’ book on fPortfolio (archived here). The bug was described nearly 10 years ago and it still hasn’t been fixed, so I figured we could turn to the community to help address it. My goal in creating this issue is to summarize what is known, identify potential issues with the package code, and hopefully get help to address it.
Here are the StackOverflow issues on this issue:
- Setting target risk in R package fPortfolio (archived here)
- Using target risk or target return in R package fPortfolio (archived here)
Here is a minimal reproducible example:
library("quantmod")
#> Warning: package 'quantmod' was built under R version 4.3.3
#> 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
#> Warning: package 'TTR' was built under R version 4.3.3
#> Registered S3 method overwritten by 'quantmod':
#> method from
#> as.zoo.data.frame zoo
library("fPortfolio")
#> Warning: package 'fPortfolio' was built under R version 4.3.3
#> Loading required package: timeDate
#> Warning: package 'timeDate' was built under R version 4.3.2
#> Loading required package: timeSeries
#> Warning: package 'timeSeries' was built under R version 4.3.3
#>
#> Attaching package: 'timeSeries'
#> The following object is masked from 'package:zoo':
#>
#> time<-
#> The following objects are masked from 'package:graphics':
#>
#> lines, points
#> Loading required package: fBasics
#> Warning: package 'fBasics' was built under R version 4.3.3
#>
#> Attaching package: 'fBasics'
#> The following object is masked from 'package:TTR':
#>
#> volatility
#> Loading required package: fAssets
#> Warning: package 'fAssets' was built under R version 4.3.3
library("dplyr")
#> Warning: package 'dplyr' was built under R version 4.3.2
#>
#> ######################### Warning from 'xts' package ##########################
#> # #
#> # The dplyr lag() function breaks how base R's lag() function is supposed to #
#> # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
#> # source() into this session won't work correctly. #
#> # #
#> # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
#> # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
#> # dplyr from breaking base R's lag() function. #
#> # #
#> # Code in packages is not affected. It's protected by R's namespace mechanism #
#> # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
#> # #
#> ###############################################################################
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:timeSeries':
#>
#> filter, lag
#> The following objects are masked from 'package:xts':
#>
#> first, last
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
# Download historical stock prices
symbols <- c("AAPL", "MSFT", "GOOGL", "AMZN", "JNJ", "JPM", "V", "PG", "XOM", "TSLA")
quantmod::getSymbols(symbols)
#> [1] "AAPL" "MSFT" "GOOGL" "AMZN" "JNJ" "JPM" "V" "PG" "XOM"
#> [10] "TSLA"
# Calculate stock returns
prices <- do.call(merge, lapply(symbols, function(sym) quantmod::Cl(get(sym))))
returns <- na.omit(TTR::ROC(prices, type = "discrete"))
returns_ts <- timeSeries::as.timeSeries(returns)
# Create portfolio
portfolioSpec <- fPortfolio::portfolioSpec()
# Identify tangency portfolio
tangencyPortfolio <- fPortfolio::tangencyPortfolio(
data = returns_ts,
spec = portfolioSpec)
# Extract optimal weights for tangency portfolio
fPortfolio::getWeights(tangencyPortfolio)
#> AAPL.Close MSFT.Close GOOGL.Close AMZN.Close JNJ.Close JPM.Close
#> 0.21136483 0.12495822 0.01223435 0.18743248 0.00000000 0.00000000
#> V.Close PG.Close XOM.Close TSLA.Close
#> 0.23495738 0.06078310 0.00000000 0.16826965
# Define target risk levels
targetRisks <- seq(0, 0.3, by = 0.01)
# Initialize storage for optimal portfolios
optimalPortfolios <- list()
optimalWeights_list <- list()
# Find optimal weightings for each target risk level
for (risk in targetRisks) {
# Create a portfolio optimization specification with the target risk
portfolioSpec <- fPortfolio::portfolioSpec()
fPortfolio::setTargetRisk(portfolioSpec) <- risk
# Solve for the maximum return at this target risk
optimal_portfolio <- fPortfolio::maxreturnPortfolio(
returns_ts,
spec = portfolioSpec)
# Store the optimal portfolio
optimalPortfolios[[as.character(risk)]] <- optimal_portfolio
# Store the optimal portfolio weights with risk level
optimal_weights <- fPortfolio::getWeights(optimal_portfolio)
optimalWeights_list[[as.character(risk)]] <- c(RiskLevel = risk, optimal_weights)
}
optimalWeightsByRisk <- dplyr::bind_rows(optimalWeights_list)
optimalWeightsByRisk
#> # A tibble: 31 × 11
#> RiskLevel AAPL.Close MSFT.Close GOOGL.Close AMZN.Close JNJ.Close JPM.Close
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 0 0.0276 0 0.0282 0.0466 0.393 0
#> 2 0.01 0.0276 0 0.0282 0.0466 0.393 0
#> 3 0.02 0.0276 0 0.0282 0.0466 0.393 0
#> 4 0.03 0.0276 0 0.0282 0.0466 0.393 0
#> 5 0.04 0.0276 0 0.0282 0.0466 0.393 0
#> 6 0.05 0.0276 0 0.0282 0.0466 0.393 0
#> 7 0.06 0.0276 0 0.0282 0.0466 0.393 0
#> 8 0.07 0.0276 0 0.0282 0.0466 0.393 0
#> 9 0.08 0.0276 0 0.0282 0.0466 0.393 0
#> 10 0.09 0.0276 0 0.0282 0.0466 0.393 0
#> # ℹ 21 more rows
#> # ℹ 4 more variables: V.Close <dbl>, PG.Close <dbl>, XOM.Close <dbl>,
#> # TSLA.Close <dbl>
sessionInfo()
#> R version 4.3.1 (2023-06-16 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 11 x64 (build 22631)
#>
#> Matrix products: default
#>
#>
#> locale:
#> [1] LC_COLLATE=English_United States.utf8
#> [2] LC_CTYPE=English_United States.utf8
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United States.utf8
#>
#> time zone: America/Chicago
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] dplyr_1.1.4 fPortfolio_4023.84 fAssets_4023.85
#> [4] fBasics_4032.96 timeSeries_4032.109 timeDate_4032.109
#> [7] quantmod_0.4.26 TTR_0.24.4 xts_0.14.0
#> [10] zoo_1.8-12
#>
#> loaded via a namespace (and not attached):
#> [1] utf8_1.2.4 generics_0.1.3 bitops_1.0-7
#> [4] robustbase_0.99-3 slam_0.1-51 lattice_0.22-6
#> [7] digest_0.6.36 magrittr_2.0.3 rneos_0.4-0
#> [10] evaluate_0.24.0 grid_4.3.1 fCopulae_4022.85
#> [13] mvtnorm_1.2-5 fastmap_1.2.0 jsonlite_1.8.8
#> [16] energy_1.7-11 fansi_1.0.6 kernlab_0.9-32
#> [19] truncnorm_1.0-9 XML_3.99-0.17 numDeriv_2016.8-1.1
#> [22] ecodist_2.1.3 mnormt_2.1.1 cli_3.6.3
#> [25] rlang_1.1.4 Rglpk_0.6-5.1 gsl_2.1-8
#> [28] reprex_2.1.1 withr_3.0.0 yaml_2.3.9
#> [31] tools_4.3.1 parallel_4.3.1 mvnormtest_0.1-9-3
#> [34] boot_1.3-30 curl_5.2.1 R6_2.5.1
#> [37] vctrs_0.6.5 stats4_4.3.1 lifecycle_1.0.4
#> [40] fs_1.6.4 MASS_7.3-60.0.1 Rsolnp_1.16
#> [43] pkgconfig_2.0.3 pillar_1.9.0 fMultivar_4031.84
#> [46] glue_1.7.0 Rcpp_1.0.13 tidyselect_1.2.1
#> [49] tibble_3.2.1 DEoptimR_1.1-3 xfun_0.46
#> [52] rstudioapi_0.16.0 knitr_1.48 spatial_7.3-17
#> [55] htmltools_0.5.8.1 igraph_2.0.3 rmarkdown_2.27
#> [58] compiler_4.3.1 quadprog_1.5-8 sn_2.1.1
#> [61] RCurl_1.98-1.16
Created on 2024-07-23 with reprex v2.1.1
Here is the package code for the maxreturnPortfolio() function:
> fPortfolio::maxreturnPortfolio
function (data, spec = portfolioSpec(), constraints = "LongOnly")
{
data = portfolioData(data, spec)
if (is.null(getTargetRisk(spec))) {
stop("Missing target risk for maximum return optimization.")
}
else {
Solver = match.fun(getSolver(spec))
portfolio = Solver(data, spec, constraints)
setWeights(spec) = portfolio$weights
setStatus(spec) = portfolio$status
Title = "Return Maximized Efficient Portfolio"
}
portfolio <- feasiblePortfolio(data, spec, constraints)
portfolio@call <- match.call()
portfolio@title <- Title
portfolio
}
It appears to call the feasiblePortfolio() function:
> fPortfolio::feasiblePortfolio
function (data, spec = portfolioSpec(), constraints = "LongOnly")
{
Data <- portfolioData(data, spec)
if (inherits(data, "fPFOLIODATA"))
data <- getSeries(Data)
assetsNames <- getUnits(Data)
Spec <- spec
Constraints <- portfolioConstraints(Data, spec, constraints)
if (is.null(getWeights(spec))) {
stop("Missing weights")
}
weights <- as.vector(getWeights(spec))
names(weights) <- assetsNames
if (inherits(getSeries(Data), "timeSeries")) {
targetReturn <- c(mean = (Data@statistics$mean %*% weights)[[1]],
mu = (Data@statistics$mu %*% weights)[[1]])
setTargetReturn(spec) <- targetReturn
Cov <- Data@statistics$Cov
cov <- sqrt((weights %*% Cov %*% weights)[[1]])
if (getType(spec) == "SPS") {
myCheck <- TRUE
funSigma <- match.fun(getObjective(spec)[1])
rcov <- funSigma(as.vector(weights))
}
else {
Sigma <- Data@statistics$Sigma
rcov <- sqrt((weights %*% Sigma %*% weights)[[1]])
}
alpha <- getAlpha(spec)
returns <- getDataPart(getSeries(Data)) %*% weights
VaR <- quantile(returns, alpha, type = 1)
CVaR <- VaR - 0.5 * mean(((VaR - returns) + abs(VaR -
returns)))/alpha
targetRisk <- c(cov, rcov, -CVaR, -VaR)
names(targetRisk) <- c("Cov", "Sigma", "CVaR", "VaR")
alpha <- getAlpha(Spec)
}
else if (inherits(getSeries(Data), "logical")) {
targetReturn <- c(mean = (Data@statistics$mean %*% weights)[[1]],
mu = NA)
setTargetReturn(spec) <- targetReturn
Cov <- Data@statistics$Cov
cov <- sqrt((weights %*% Cov %*% weights)[[1]])
targetRisk <- c(cov, NA, NA, NA)
names(targetRisk) <- c("Cov", "Sigma", "CVaR", "VaR")
alpha <- NA
}
covRiskBudgets <- (weights * Cov %*% weights)[, 1]/cov^2
names(covRiskBudgets) <- assetsNames
Portfolio <- new("fPFOLIOVAL", portfolio = list(weights = weights,
covRiskBudgets = covRiskBudgets, targetReturn = targetReturn,
targetRisk = targetRisk, targetAlpha = alpha, status = getStatus(spec)))
new("fPORTFOLIO", call = match.call(), data = Data, spec = Spec,
constraints = Constraints, portfolio = Portfolio, title = "Feasible Portfolio",
description = description())
}
As far as I can tell based on the code (and I’m no expert), the feasiblePortfolio() function seems to create and assign the target risk, not to accept it as an input from the user (i.e., targetRisk <- c(cov, rcov, -CVaR, -VaR)), although I’m not sure. As one SO user noted, “[It] seems to have something to do with the efficientPortfolio() function and the solver used in spec. Apparently no matter what targetRisk you input, the solver inputs an objective that corresponds to the Covyou get when you print the portfolio”.
Any help in helping to identify and fix this issue would be greatly appreciated.