@Axel Vogt I thought you would never ask (not marriage) !

" Edited: I do not understand why the point ~ (1,1) is 'optimal' for you and
what is the meaning of the falling line"

We are interested in finding the minimum value of the portfolio y[1]+y[2]
on the 0.05 probability multivariate Quantile ie what is the maxium portfolio
loss on the 5% level. That point seams to be approx [-1,-1]. Note that:

i) I have not yet confirmed if y[1]+y[2] is constant on the 0.05 prob Quantile.
ii) if the optimal point indeed is y[1]=y[2]

Any suggestions on how we can find such point or any otherway of solving
the problem most appreciated

This is too complicated. If I instead of trying to find the 5% Quantile
of two multivariate random variables try to find
the expected value of the sum of two multivariate random variables (x+y)
then I indeed can find a solution.
 

restart:
with(Optimization):
with(Statistics):
with(plots):

rho := 0:
mu[x] := 0:
mu[y] := 0:
sigma[x] := 1:
sigma[y] := 1:

pdf := exp(-((x-mu[x])^2/sigma[x]^2-2*rho*(x-mu[x])*(y-mu[y])/(sigma[x]*sigma[y])+(y-mu[y])^2/sigma[y]^2)/(2*(1-rho^2)))/(2*Pi*sigma[x]*sigma[y]*sqrt(1-rho^2));

Maximize(pdf*(x+y), x = -4 .. 4, y = -4 .. 4);
Minimize(pdf*(x+y), x = -4 .. 4, y = -4 .. 4);
plot3d(pdf*(x+y), x = -4 .. 4, y = -4 .. 4, axes = boxed);

[0.136517362297204348, [x = 0.707106781209232316, y = 0.707106781209232316]]
[-0.136517362297204348, [x = -0.707106781209232316, y = -0.707106781209232316]]

This is too complicated. If I instead of trying to find the 5% Quantile
of two multivariate random variables try to find
the expected value of the sum of two multivariate random variables (x+y)
then I indeed can find a solution.
 

restart:
with(Optimization):
with(Statistics):
with(plots):

rho := 0:
mu[x] := 0:
mu[y] := 0:
sigma[x] := 1:
sigma[y] := 1:

pdf := exp(-((x-mu[x])^2/sigma[x]^2-2*rho*(x-mu[x])*(y-mu[y])/(sigma[x]*sigma[y])+(y-mu[y])^2/sigma[y]^2)/(2*(1-rho^2)))/(2*Pi*sigma[x]*sigma[y]*sqrt(1-rho^2));

Maximize(pdf*(x+y), x = -4 .. 4, y = -4 .. 4);
Minimize(pdf*(x+y), x = -4 .. 4, y = -4 .. 4);
plot3d(pdf*(x+y), x = -4 .. 4, y = -4 .. 4, axes = boxed);

[0.136517362297204348, [x = 0.707106781209232316, y = 0.707106781209232316]]
[-0.136517362297204348, [x = -0.707106781209232316, y = -0.707106781209232316]]

Not really but thanx anyway :-)

Not really but thanx anyway :-)

Joe, yes you are right. I forgot with(Statistics) in the begining.
However, I still get an error message (see below). W is a vector column ie <w[1],w[2],w[3]..w[n]>
which means that data1 is a vector column of the form:

R[1,1]*w[1]+R[1,2]*w[2]+.....R[1,n]*w[n]    # Portfolio Return t=1
R[2,1]*w[1]+R[2,2]*w[2]+.....R[2,n]*w[n]    # Portfolio Return t=2

in other words data1 is a vector column containing the Portfolio Return from t=1 to t=n.
We want to find the weights w[i] that minimize the 5th percentile of such a time series.


Error, (in Statistics:-Percentile) expected 1-dimensional Array, received: Array(1..15, 1..1, {(1, 1) = 5.*w[1]+13.*w[2]+4.*w[3]-15.*w[4]-2.*w[5]+4.*w[6]+11.*w[7]+11.*w[8]+11.*w[9]+3.*w[10]+11.*w[11]+2.*w[12]+10.*w[13]+6.*w[14]-2.*w[15]+15.*w[16]+1.*w[17]+13.*w[19]+3.*w[20], (2, 1) = 11.*w[1]+11.*w[2]-6.*w[3]-11.*w[4]+4.*w[6]-3.*w[7]+2.*w[8]-1.*w[10]+6.*w[12]-3.*w[13]-11.*w[14]-9.*w[15]-12.*w[16]-3.*w[17]+10.*w[18]-11.*w[19]+11.*w[20], (3, 1) = 1.*w[1]+6.*w[2]+4.*w[3]-2.*w[4]-3.*w[5]-13.*w[6]+1.*w[7]-4.*w[8]+3.*w[9]-3.*w[10]+1.*w[11]-10.*w[12]+2.*w[13]+3.*w[14]-14.*w[15]-5.*w[16]-14.*w[17]+9.*w[18]+9.*w[19]-4.*w[20], (4, 1) = -3.*w[1]-1.*w[2...

Joe, yes you are right. I forgot with(Statistics) in the begining.
However, I still get an error message (see below). W is a vector column ie <w[1],w[2],w[3]..w[n]>
which means that data1 is a vector column of the form:

R[1,1]*w[1]+R[1,2]*w[2]+.....R[1,n]*w[n]    # Portfolio Return t=1
R[2,1]*w[1]+R[2,2]*w[2]+.....R[2,n]*w[n]    # Portfolio Return t=2

in other words data1 is a vector column containing the Portfolio Return from t=1 to t=n.
We want to find the weights w[i] that minimize the 5th percentile of such a time series.


Error, (in Statistics:-Percentile) expected 1-dimensional Array, received: Array(1..15, 1..1, {(1, 1) = 5.*w[1]+13.*w[2]+4.*w[3]-15.*w[4]-2.*w[5]+4.*w[6]+11.*w[7]+11.*w[8]+11.*w[9]+3.*w[10]+11.*w[11]+2.*w[12]+10.*w[13]+6.*w[14]-2.*w[15]+15.*w[16]+1.*w[17]+13.*w[19]+3.*w[20], (2, 1) = 11.*w[1]+11.*w[2]-6.*w[3]-11.*w[4]+4.*w[6]-3.*w[7]+2.*w[8]-1.*w[10]+6.*w[12]-3.*w[13]-11.*w[14]-9.*w[15]-12.*w[16]-3.*w[17]+10.*w[18]-11.*w[19]+11.*w[20], (3, 1) = 1.*w[1]+6.*w[2]+4.*w[3]-2.*w[4]-3.*w[5]-13.*w[6]+1.*w[7]-4.*w[8]+3.*w[9]-3.*w[10]+1.*w[11]-10.*w[12]+2.*w[13]+3.*w[14]-14.*w[15]-5.*w[16]-14.*w[17]+9.*w[18]+9.*w[19]-4.*w[20], (4, 1) = -3.*w[1]-1.*w[2...

@acer

i) "Perhaps a timing comparison is in order here. The Optimization package does these QP problems with hundreds of variables in a few seconds on a fast machine. Marcus gave one with, and one without, general constraints. I wonder how DirectSearch performs on them when passed its method=quadratic option."

The performance is not good :-(. I wish I could say it was. I have done it and it took a long long long time
with 600 variables with constraints. I would upload the code but unfortunatly I have deleted it.

ii) "Note that, as yet, Marcus has not answered posed questions about what size problems he is aiming for."

I am aiming for the sky ie 1000's of variables he he. No honestly I prefer quality over enormous datasets (but it would be cool). For example getting allocations like -1,+1,0.2,-1,+1 for 500 variables just makes you want to cry for the max risk adjusted return long-short QP matrix form algorithm. Long only sucks also.

@Pseudomodo  I ran this "2-stage method" but it is much slower than QP matrix form
so I think I agree with Acer on this one. A global solver like with(DirecSearch) might
not be optimal to use for these kind of problems.

@Markiyan Hirnyk your persistence has served you well! Exellent, thank you for sharing.

ahhhh thank you of course. Ok, last question I promise :-)
When I run the minimize portfolio variance code (below) both for algebraic and
matrix form it produces the same output (yipee) however both results in an error message

Warning, necessary conditions met but sufficient conditions not satisfied

Is there anything we can do about that?! I mean to remove it...?!

# Minimize Portfolio Variance - Algebraic and Matrix Form
# Minimize W'.Cov.W

restart:
with(LinearAlgebra):
with(Optimization):

N := 600:
U := RandomMatrix(N, outputoptions = [datatype = float[8], shape = triangular[upper]]):
Cov := Matrix(U.U^%T, shape = symmetric):

EV := RandomVector(N, outputoptions = [datatype = float[8]]):
A := Matrix([[seq(1, i = 1 .. N)]], datatype = float[8]):
b := Vector([1], datatype = float[8]):
bl := Vector(N, fill = 0, datatype = float[8]):
bu := Vector(N, fill = 1, datatype = float[8]):
sol1 := QPSolve(2*Cov, [A, b], [bl, bu]);

W := Vector(N, symbol = w):
S := Vector(N, fill = 1, datatype = float[8]):

z := QPSolve(W^%T.Cov.W, {S.W <= 1}, seq(W[i] = bl[i] .. bu[i], i = 1 .. N)):
sol2 := [z[1], convert(map(proc (x) options operator, arrow; rhs(x) end proc, z[2]), Vector)];


Warning, necessary conditions met but sufficient conditions not satisfied
                [                        [ 1 .. 600 Vector[column] ]]
                [                   -59  [ Data Type: float[8]     ]]
        sol1 := [2.19215436624030 10   , [ Storage: rectangular    ]]
                [                        [ Order: Fortran_order    ]]

Warning, necessary conditions met but sufficient conditions not satisfied
                [                        [ 1 .. 600 Vector[column] ]]
                [                   -59  [ Data Type: anything     ]]
        sol2 := [2.19215436624030 10   , [ Storage: rectangular    ]]
                [                        [ Order: Fortran_order    ]]

ahhhh thank you of course. Ok, last question I promise :-)
When I run the minimize portfolio variance code (below) both for algebraic and
matrix form it produces the same output (yipee) however both results in an error message

Warning, necessary conditions met but sufficient conditions not satisfied

Is there anything we can do about that?! I mean to remove it...?!

# Minimize Portfolio Variance - Algebraic and Matrix Form
# Minimize W'.Cov.W

restart:
with(LinearAlgebra):
with(Optimization):

N := 600:
U := RandomMatrix(N, outputoptions = [datatype = float[8], shape = triangular[upper]]):
Cov := Matrix(U.U^%T, shape = symmetric):

EV := RandomVector(N, outputoptions = [datatype = float[8]]):
A := Matrix([[seq(1, i = 1 .. N)]], datatype = float[8]):
b := Vector([1], datatype = float[8]):
bl := Vector(N, fill = 0, datatype = float[8]):
bu := Vector(N, fill = 1, datatype = float[8]):
sol1 := QPSolve(2*Cov, [A, b], [bl, bu]);

W := Vector(N, symbol = w):
S := Vector(N, fill = 1, datatype = float[8]):

z := QPSolve(W^%T.Cov.W, {S.W <= 1}, seq(W[i] = bl[i] .. bu[i], i = 1 .. N)):
sol2 := [z[1], convert(map(proc (x) options operator, arrow; rhs(x) end proc, z[2]), Vector)];


Warning, necessary conditions met but sufficient conditions not satisfied
                [                        [ 1 .. 600 Vector[column] ]]
                [                   -59  [ Data Type: float[8]     ]]
        sol1 := [2.19215436624030 10   , [ Storage: rectangular    ]]
                [                        [ Order: Fortran_order    ]]

Warning, necessary conditions met but sufficient conditions not satisfied
                [                        [ 1 .. 600 Vector[column] ]]
                [                   -59  [ Data Type: anything     ]]
        sol2 := [2.19215436624030 10   , [ Storage: rectangular    ]]
                [                        [ Order: Fortran_order    ]]

Now, I have run into a new problem. The algebraic and matrix form of the minimize portfolio variance
for a given expected return produces different results. Which indicate to me that something is wrong.
What is wrong this time?!

restart:
with(LinearAlgebra):
with(Optimization):

N := 600:
U := RandomMatrix(N, outputoptions = [datatype = float[8], shape = triangular[upper]]):
Cov := Matrix(U.U^%T, shape = symmetric):

EV := RandomVector(N, outputoptions = [datatype = float[8]]):
pr := Vector([.6*max(EV)], datatype = float[8]):
A := Matrix([[seq(1, i = 1 .. N)]], datatype = float[8]):
b := Vector([1], datatype = float[8]):
bl := Vector(N, fill = 0, datatype = float[8]):
bu := Vector(N, fill = 1, datatype = float[8]):

EV2 := Transpose(Matrix(EV, datatype = float[8])):

QPSolve(Cov, [A, b, EV2, pr], [bl, bu]);

W := Vector(N, symbol = w):
S := Vector(N, fill = 1, datatype = float[8]):

z:=QPSolve(W^%T.Cov.W, {EV.W = pr[1], S.W <= 1}, seq(W[i] = bl[i] .. bu[i], i = 1 .. N)):

[z[1], convert(map(proc (x) options operator, arrow; rhs(x) end proc, z[2]), Vector)];


            
               [                                          [ 1 .. 600 Vector[column] ]]
               [                                          [ Data Type: float[8]     ]]
               [37.9169642447005,         [ Storage: rectangular    ]]
               [                                         [ Order: Fortran_order    ]]
               
               [                                         [ 1 .. 600 Vector[column] ]]
               [                                         [ Data Type: anything     ]]
               [75.8339284894005,          [ Storage: rectangular    ]]
               [                                         [ Order: Fortran_order    ]]

Now, I have run into a new problem. The algebraic and matrix form of the minimize portfolio variance
for a given expected return produces different results. Which indicate to me that something is wrong.
What is wrong this time?!

restart:
with(LinearAlgebra):
with(Optimization):

N := 600:
U := RandomMatrix(N, outputoptions = [datatype = float[8], shape = triangular[upper]]):
Cov := Matrix(U.U^%T, shape = symmetric):

EV := RandomVector(N, outputoptions = [datatype = float[8]]):
pr := Vector([.6*max(EV)], datatype = float[8]):
A := Matrix([[seq(1, i = 1 .. N)]], datatype = float[8]):
b := Vector([1], datatype = float[8]):
bl := Vector(N, fill = 0, datatype = float[8]):
bu := Vector(N, fill = 1, datatype = float[8]):

EV2 := Transpose(Matrix(EV, datatype = float[8])):

QPSolve(Cov, [A, b, EV2, pr], [bl, bu]);

W := Vector(N, symbol = w):
S := Vector(N, fill = 1, datatype = float[8]):

z:=QPSolve(W^%T.Cov.W, {EV.W = pr[1], S.W <= 1}, seq(W[i] = bl[i] .. bu[i], i = 1 .. N)):

[z[1], convert(map(proc (x) options operator, arrow; rhs(x) end proc, z[2]), Vector)];


            
               [                                          [ 1 .. 600 Vector[column] ]]
               [                                          [ Data Type: float[8]     ]]
               [37.9169642447005,         [ Storage: rectangular    ]]
               [                                         [ Order: Fortran_order    ]]
               
               [                                         [ 1 .. 600 Vector[column] ]]
               [                                         [ Data Type: anything     ]]
               [75.8339284894005,          [ Storage: rectangular    ]]
               [                                         [ Order: Fortran_order    ]]

That is an excellent answer on so many levels! Thank you for that.

i) Both algorithms (algebraic and matrix form) produces the same output
and the same allocations. Check one.

ii) The matrix QP form runs faster. So matrix form is better to use. Check two.

iii) The [-EV,2*Cov] did not make any sense when I looked at it first but when I saw that the
default form for QP was  c' x + (1/2) x' H x it all made sense. Which means  x' COV x - Ev'x which is exactly
the form we want. Again thank you for your expertise :-)