restart:
with(LinearAlgebra):
with(ArrayTools):
with(Statistics):
randomize():
with(plots):
with(combinat):

n := 100:
nstock := 7:
corr := .8:

R := Matrix(nstock, nstock, proc (i, j) options operator, arrow; `if`(i = j, 1, corr) end proc):
CD := Matrix(LUDecomposition(evalf(R), 'method' = 'Cholesky'), datatype = float[8]):

ev := [seq((1/5)*(rand(-3 .. 4))(), i = 1 .. nstock)]:
st := [seq((rand(1 .. 2))(), i = 1 .. nstock)]:
A := Matrix([seq([Sample(RandomVariable(Normal(ev[i], st[i])), n)], i = 1 .. nstock)]):
B := Matrix(Transpose(Multiply(CD, A)));

ap := choose([seq(i, i = 1 .. nstock), seq(-i, i = 1 .. nstock)])[2 .. numbcomb([seq(i, i = 1 .. nstock), seq(-i, i = 1 .. nstock)])]:

for i to nops(ap) do

PR[i] := add(`if`(ap[i][j] > 0, 1, -1)*Array(Column(B, ap[i][j])), j = 1 .. nops(ap[i])):

EVV[i] := ExpectedValue(PR[i]):
STT[i] := StandardDeviation(PR[i])

end do:

pointplot([seq([STT[i], EVV[i]], i = 1 .. nops(ap))], color = COLOR(RGB, rand()/10^12, rand()/10^12, rand()/10^12), thickness = 3, symbol = solidcircle, symbolsize = 10);