The parma package offers a set of models and methods for optimal portfolio allocation based on scenario optimization. It has certain unique features over competing products, such as the ability to optimize directly the risk to reward ratio using fractional programming, and offers LP, MILP, QP, NLP, and GNLP formulations for different problems. For convex NLP problems, all functions and constraints are calculated making use of analytic derivatives for confident optimization, and for certain NLP problems which were previously hard to optimize because of discontinuous functions in their makeup, smooth approximations are used which offer very high accuracy.
The parmaspec function
The first step in optimizing is defining a specification of the problem using parmaspec:
library(parma, quietly = TRUE, verbose = FALSE) args(parmaspec) ## function (scenario = NULL, probability = NULL, S = NULL, benchmark = NULL, ## benchmarkS = NULL, forecast = NULL, target = NULL, targetType = c('inequality', ## 'equality'), risk = c('MAD', 'MiniMax', 'CVaR', 'CDaR', ## 'EV', 'LPM', 'LPMUPM'), riskType = c('minrisk', 'optimal'), ## options = list(alpha = 0.05, threshold = 999, moment = 1, ## lmoment = 1, umoment = 1, lthreshold = -0.01, uthreshold = 0.01), ## LB = NULL, UB = NULL, budget = 1, leverage = NULL, ineqfun = NULL, ## ineqgrad = NULL, eqfun = NULL, eqgrad = NULL, uservars = list(), ## ineq.mat = NULL, ineq.LB = NULL, ineq.UB = NULL, eq.mat = NULL, ## eqB = NULL, max.pos = NULL, asset.names = NULL, ...) ## NULL
There is a choice of either providing a scenario matrix or a covariance matrix S, in which case the risk will automatically be set to “EV” (Variance). There is a choice of 6 deviation/risk measures (the “LPMUPM” is not a separate risk measure but one with a reward measure -UPM- different from the typical linear forecast vector), 2 methods for optimizing (minimum risk or the fractional optimal risk to reward), and 2 methods for defining the reward constraint (as an inequality or equality). A benchmark vector may be provided in which case benchmark relative optimization may be used, for which it is usual to set the lower (LB) and upper (UB) bounds to the allowable positive and negative deviations from the benchmark, with a budget constraint of zero so that the resulting optimal weights are ‘active’ weights. Similarly for the target and forecast, these should then be considered active, and appropriate values passed to them. The remaining options define constraints, the type of which may limit the problem to a specific form. For example, many of the measures may be solved using either LP or NLP methods, but the moment leverage is used (for example) with short selling, the problem is automatically defined as NLP and a smooth approximation to the absolute value function used together with its derivative (thus avoiding the trimabillity assumptions and doubling of the problem space in Jacobs, Levy and Markowitz (2006)). In a separate post the use of custom constraints for NLP is demonstrated, examples of which are already available in the parma.tests folder of the package.
The following example demonstrates how to define and solve a typical problem. The use of the static data matrix is only for illustration purposes and is strongly discouraged in general. The rmgarch package now contains a set of functions called fscenario and fmoments for creating forecast scenarios and moments for use with parma.
library(quantmod) data(etfdata) # optionally set a timezone to avoid getting times in the index: # Sys.setenv(TZ = "GMT") R = ROC(etfdata, na.pad=FALSE) spec = parmaspec(scenario = R, forecast = colMeans(R), risk = 'MAD', target = mean(colMeans(R)), targetType = 'equality', riskType = 'minrisk', LB = rep(0, 15), UB = rep(1, 15), budget = 1) show(spec) ## ## +---------------------------------+ ## | PARMA Specification | ## +---------------------------------+ ## No.Assets : 15 ## Problem : LP,NLP,GNLP ## Input : Scenario ## Risk Measure : MAD ## Objective : minrisk
The problem can be solved as either an LP, NLP or GNLP problem (the latter should not be used unless the problem is LPMUPM or any other which has custom constraints which are known to be non-convex).
The parmasolve function
This takes as arguments the portfolio specification, of class parmaSpec, the type of problem formulation (if a choice exists), an optional solver control list (for the nloptr solver), and some optional parameters for use if the problem type is GNLP.
sol_lp = parmasolve(spec, type = 'LP') sol_nlp = parmasolve(spec, type = 'NLP')
The resulting solution objects, of class parmaPort, have certain methods, including a summary, to extract their values and displayed below:
out = cbind(weights(sol_lp), weights(sol_nlp)) out = rbind(out, unname(cbind(parmarisk(sol_lp), parmarisk(sol_nlp)))) out = rbind(out, cbind(parmareward(sol_lp), parmareward(sol_nlp))) out = rbind(out, cbind(as.numeric(tictoc(sol_lp)), as.numeric(tictoc(sol_nlp)))) rownames(out)[16:18] = c('risk', 'reward', 'elapsed(secs)') colnames(out) = c('MAD[LP]', 'MAD[NLP]') print(round(out, 5)) ## MAD[LP] MAD[NLP] ## IWF 0.20660 0.20660 ## IWD 0.00000 0.00000 ## IWO 0.00000 0.00000 ## IWN 0.00000 0.00000 ## EEM 0.00000 0.00000 ## TLT 0.61768 0.61768 ## EWC 0.09973 0.09973 ## EWA 0.00000 0.00000 ## EWJ 0.00000 0.00000 ## EWG 0.00000 0.00000 ## EWL 0.01130 0.01130 ## EWQ 0.00000 0.00000 ## EWU 0.00000 0.00000 ## EPP 0.04260 0.04260 ## EZA 0.02210 0.02210 ## risk 0.00434 0.00434 ## reward 0.00032 0.00032 ## elapsed(secs) 2.20215 1.05469
The parmafrontier function
This example demonstrates the solution to the fractional programming problem and the use of the parmafrontier function.
# create a cluster object to speed up the frontier calculations: library(parallel) cl = makePSOCKcluster(10) f = colMeans(R) spec1 = parmaspec(scenario = R, target = 0, targetType = 'equality', forecast = f, risk = 'CVaR', riskType = 'minrisk', options = list(alpha = 0.05), LB = rep(0,15), UB = rep(1, 15), budget = 1) flp = parmafrontier(spec1, n.points = 100, type = 'LP', cluster = cl) r1 = flp[, 'CVaR'] rw1 = flp[, 'reward'] fidx = which(r1/rw1 == min(r1/rw1)) # Also calculate the NLP representation of the problem fnlp = parmafrontier(spec1, n.points = 100, type = 'NLP', cluster = cl) r2 = fnlp[, 'CVaR'] rw2 = fnlp[, 'reward'] # Now calculate instead directly the optimal spec2 = parmaspec(scenario = R, forecast = f, risk = 'CVaR', riskType = 'optimal', options = list(alpha = 0.05), LB = rep(0, 15), UB = rep(1, 15), budget = 1) opt_lp = parmasolve(spec2, type = 'LP')
The plot below demonstrates the high accuracy of the NLP approximation of the CVaR problem, and the use of fractional programming to obtain the optimal value rather than having to ‘find’ it on the frontier. More details on fractional programming may be found in Charnes and Cooper (1962), and more recently in Stoyanov, Rachev and Fabozzi (2007).
plot(r1, rw1, type = 'p', col = 'steelblue', xlab = 'E[Risk]', ylab = 'E[Return]', main = 'CVaR Frontier') lines(r2, rw2, col = 'tomato1') points(r1[fidx], rw1[fidx], col = 'orange', pch = 4, lwd = 3) points(parmarisk(opt_lp), parmareward(opt_lp), col = 'yellow', pch = 6, lwd = 3) legend('topleft', c('LP Frontier', 'NLP Frontier', paste('Frontier Optimal (', round(r1[fidx]/rw1[fidx], 3), ')', sep = ''), paste('Fractional Optimal (', round(parmarisk(opt_lp)/parmareward(opt_lp), 3), ')', sep = '')), col = c('steelblue', 'tomato1', 'orange', 'yellow'), pch = c(1, -1, 4, 6), lwd = c(1, 1, 3, 3), lty = c(0, 1, 0, 0), bty = 'n', cex = 0.9) grid() mtext(paste('parma package', sep = ''), side = 4, adj = 0, padj = -0.6, col = 'grey50', cex = 0.7)
Charnes, A., & Cooper, W. (1962). Programming with linear fractional functionals. Nav. Res. Logist. Q., 9, 181-186.
Jacobs, B. I., Levy, K. N., & Markowitz, H. M. (2006). Trimability and fast optimization of long-short portfolios. Financial Analysts Journal 62(2), 36-46.
Stoyanov, S. V., Rachev, S. T., & Fabozzi, F. J. (2007). Optimal Financial Portfolios. Applied Mathematical Finance, 14(5), 401-436.