risk parity portfolio from scenario optimization

Expanding my comment on non-standard risk measures: under risk parity/ERC, you look at the marginal contributions to risk, which means the first derivatives of a function $f$ that maps weights (and data) into a single number, the risk. If $f$ is the standard deviation, it can be conveniently and efficiently computed from the covariance matrix. But after all, you look for first derivatives, which you can approximate by finite differences. The downside, of course, is a higher computational cost, but there might be no other alternative way to compute those contributions, and you gain flexibility. Specifically, if $R$ is the returns matrix, and $w$ are the weights, then the $i$-th weight's derivative can be approximated as $(f(R(w^*_i)) - f(Rw))/ h$. The symbol $h$ is a small number: for a forward difference, as here, typically of the order of $10^{-8}$ or similar; when in doubt, use a larger number. $w^*_i$ is the weights vector with the $i$-th weight increased by $h$.

Computational 'proof': I use R to simulate 20 return scenarios for 5 assets.

library("FRAPO")  ## for function 'mrc'
library("NMOF")   ## for function 'randomReturns'

na <- 5
R <- randomReturns(na = na, ns = 20, sd = 0.02, rho = 0.5)
##          [,1]    [,2]    [,3]    [,4]   [,5]
##  [1,]  0.0381  0.0320  0.0126  0.0161  0.015
##  [2,]  0.0117  0.0028  0.0367 -0.0100 -0.004
##  [3,] -0.0258 -0.0006 -0.0124  0.0017 -0.005
##  [4,]  0.0309  0.0182  0.0064  0.0197  0.014
##  [5,] -0.0202 -0.0259 -0.0457  0.0064 -0.023
##  [6,] -0.0036  0.0018 -0.0345 -0.0149 -0.019
## ....


## equal-weight portfolio
w <- rep(1/na, na)


## risk function: input is a vector of _portfolio_
## returns under the scenarios
risk <- function(Rw, ...) {
    sqrt(var(c(Rw)))
}

risk(R %*% w)
## [1] 0.01432615   ## vol of the equal-weight portfolio

The function mrc_fd implements the above formula for a finite difference:

## approximate first-derivatives of risk with respect
## to weights via forward difference
mrc_fd <- function(w, R, risk, ..., h = 1e-8) {
    Rw <- R %*% w
    dR <- R*h

    ans <- numeric(length(w))
    for (i in seq_along(w)) {    
        ans[i] <- w[i] * ( risk(Rw + dR[, i]) - risk(Rw) ) / h
    }
    ans
}

Compare the results with the implementation of marginal risk contribution in package FRAPO:

FRAPO::mrc(w, var(R), FALSE)
## [1] 0.00301 0.00200 0.00349 0.00240 0.00343

mrc_fd(w = w, R = R, risk = risk)
## [1] 0.00301 0.00200 0.00349 0.00240 0.00343

The numbers are the same. But where mrc computes them off the covariance matrix, function mrc_fd directly works with the scenarios, and so would work with other definitions of risk as well.

Concerning methods: I like local-search methods, which (like scenario optimization) are extremely flexible. Essentially, these method start with random solution, and then iteratively change this solution by increasing/decreasing weights in small steps, until they have arrived at a good solution (or some computational limit is reached, such as a fixed number of iterations). Threshold Accepting (TA), for instance, is such a method, and it is also implemented in the NMOF package. TA is also an early, perhaps the earliest derivative-free method that was used for portfolio selection, as described in

@ARTICLE{Dueck1992,
  author       = {Gunter Dueck and Peter Winker},
  title        = {New Concepts and Algorithms for Portfolio Choice},
  journal      = {Applied Stochastic Models and Data Analysis},
  year         = 1992,
  volume       = 8,
  pages        = {159--178},
  number       = 3
}

(Disclosure: I am the maintainer of package NMOF.)