You could try a heuristic approach.
The problem can be split into two nested optimisations: i) in the inner optimisation, given a set of selected assets, compute mean--variance efficient weights; ii) in the outer optimisation, you iterate through combinations of assets.
The inner optimisation can be solved via a quadratic
programme (QP). For the outer optimisation, you could
use a local-search technique. To sketch how such a
method could be applied, I create some random data. I
will assume that for every asset, there are 60 return
observations (for intuition, think of them
as monthly returns).
set.seed(75578)
I <- cbind(0:15*1440+1, 1:16*1440)
rownames(I) <- LETTERS[1:nrow(I)]
colnames(I) <- c("first", "last")
R <- rnorm(60*max(I), sd = 0.02)
dim(R) <- c(60, max(I))
The returns matrix R is of size 60 times 23040, as per your comment: 16 countries with 1440 assets each. The
columns of R are assigned to countries A, B, ... and the
start and end indices for each country are in matrix I. This latter matrix looks like this:
first last
A 1 1440
B 1441 2880
C 2881 4320
D 4321 5760
E 5761 7200
F 7201 8640
G 8641 10080
H 10081 11520
I 11521 12960
J 12961 14400
K 14401 15840
L 15841 17280
M 17281 18720
N 18721 20160
O 20161 21600
P 21601 23040
I collect all these data in a single list, which makes
it easier to pass the information to a function.
Data <- list(I = I,
nc = nrow(I),
R = R)
Now suppose we have some set of assets:
x <- Data$I[,1]
x looks as follows:
A B C D E F G H
1 1441 2881 4321 5761 7201 8641 10081
I J K L M N O P
11521 12961 14401 15841 17281 18721 20161 21601
To compute the solution to the inner problem, I use the
function mvar. To keep things simple, it computes the long-only
minimum-variance portfolio. The function makes use of the NMOF package. (I am the author of the package. You will need the latest version, which can be installed as shown below or from GitHub; with the CRAN version, you will have to write NMOF:::minvar).
## install.packages('NMOF', type = 'source',
## repos = c('http://enricoschumann.net/R',
## getOption('repos')))
require("NMOF")
mvar <- function(x, Data) {
cv <- cov(Data$R[ ,x])
w <- minvar(cv, wmin = 0, wmax = 0.5)
c(sqrt(w %*% cv %*% w))
}
We check how good the portfolio x is by computing its objective function value.
mvar(x, Data)
which, with the given seed, is
0.003645508
Actually, before we use a local search, let us try
a 'constructive' solution as a benchmark: select the asset with the
lowest volatility from every country, and then compute the
minimum-variance weights for this selection.
x_sort <- numeric(Data$nc)
for (i in 1:nrow(Data$I)) {
cols <- Data$I[i,1]:Data$I[i,2]
x_sort[i] <- cols[order(apply(Data$R[ ,cols ], 2, sd))[1]]
}
We can compute this solution's objective function value
mvar(x_sort, Data)
which is better:
0.002234832
Note that the volatility (0.2%) is an
order of magnitude lower than the average volatility in the
random data set, which was 2%.
Now, in a local search, we start with some solution and
try to improve it iteratively. For this, we need a
neighbourhood function that takes a solutions and
returns a slightly changed solution. Here is one way to
write such a function.
nb <- function(x, Data) {
## randomly pick one country
i <- sample(Data$nc, 1)
## randomly pick one asset
x[i] <- sample(Data$I[i,1]:Data$I[i,2], 1)
x
}
Now I try two local-search algorithms: a simple
stochastic local search (LSopt) and Threshold Accepting (TAopt). I
use the implementations in the NMOF package. The difference between the two algorithms is that LSopt will
only accept neighbours that are not worse than the
current solution, whereas TAopt is more forgiving and will
even accept solutions that are slightly worse than the
current solution. In this way, TAopt may escape from local
minima. Both methods use 100000 iterations.
steps <- 100000
sol_LS <- LSopt(mvar, list(x0 = Data$I[,1],
neighbour = nb,
nS = steps),
Data = Data)
sol_TA <- TAopt(mvar, list(x0 = Data$I[,1],
neighbour = nb,
nS = steps/10),
Data = Data)
We can compare the results of both methods with our
benchmark solution, which I add as a blue horizontal
line.
par(mar = c(5,5,1,1), las = 1, mgp = c(3,0.25,0), tck = 0.01)
plot(x = seq(1, steps, length.out = 1000),
y = sol_TA$Fmat[seq(1, steps, length.out = 1000),2], log = "y",
xlab = "iteration", ylab = "portfolio vol", type = "l")
lines(x = seq(1, steps, length.out = 1000),
y = sol_LS$Fmat[seq(1, steps, length.out = 1000),2], col = grey(.4))
abline(h = mvar(x_sort, Data), col = "blue")
legend("topright",
legend = c("Local Search", "Threshold Accepting"),
col = c(grey(.4), "black"), lwd = 4, lty =1)
text(x = steps*0.8, y = mvar(x_sort, Data), pos = 3,
"lowest-vol asset per country", col = "blue")
LSopt (grey) and TAopt (black) both achieve solutions of around 0.1%, so substantially better than the constructive solution. Threshold Accepting performs slightly better than a simple local search.
