ok, thanx for your input ! sounds complicated :-)

I would have hoped that I could simply click on a button. Voila !

Here is your large data set mdb file converted to Maple Mr researcher.......... 

I have no idea...... which is best and easies?

maybe I am a bit of a drama queen... it turned out not to be that hard after you got over the mental block of having to test all

these different permutations...... I think the below code will do the job

#############################

with(combinat):
P := 3:  # Portfolio Size
A := permute(nstock, P):
g := nops(permute(nstock, P));
for i to g do
rr := << Column(r,A[i,1]) >|< Column(r,A[i,2]) >|< Column(r,A[i,3]) >>:
S[i] := ExpectedValue(Column(r, A[i, 1])+Column(r, A[i, 2])+Column(r, A[i, 3]))/sqrt(add(i, `in`(i, CovarianceMatrix(rr))))
end do:

plot([seq(i, i = 1 .. g)], [seq(S[i], i = 1 .. g)], axis = [gridlines = [30, color = blue]], labels = ["Permutations", "Portfolio Sharpe"], font = [Times, Roman, 14]);
 [seq(S[i], i = 1 .. g)];

#############################

But I still would like to have some comments... I am not sure that everything is correct....

I have made some minor changes to the code. The more I play around with it the more I realize that I struck a gold mine.

I am quite impressed by the speed of the algorithm and Maple, 1300 permutations in a couple of seconds

########################################################

with(combinat):
P := 3:
A := permute(nstock, P):
g := nops(permute(nstock, P)):

for i to g do
rr := <<Column(r,A[i,1])> |  <Column(r,A[i,2])> | <Column(r,A[i,3])> >:
S[i] := ExpectedValue(Column(r, A[i, 1])+Column(r, A[i, 2])+Column(r, A[i, 3]))/sqrt(add(i, `in`(i, CovarianceMatrix(rr))))

end do:

plot([seq(i, i = 1 .. g)], [seq(S[i], i = 1 .. g)], axis = [gridlines = [30, color = blue]], labels = ["Permutations", "Portfolio Sharpe"], font = [Times, Roman, 14]);
"Max Portfolio Sharpe" = max(seq(S[i], i = 1 .. g));

for i to g do
if S[i] = max(seq(S[i], i = 1 .. g)) then W[i] := i else W[i] := 0
end if end do:
WW := RemoveInRange([seq(W[i], i = 1 .. g)], 0 .. 1):

for i to nops(WW) do
QQ[i] := A[WW[i]]
end do:
QQQ := [seq(QQ[i], i = 1 .. nops(WW))][1];


RB := <<Column(r,QQQ[1])> |  <Column(r,QQQ[2])> | <Column(r,QQQ[3])> >:

"Sharpe Ratio Portfolio" = ExpectedValue([seq(r[i, QQQ[1]]+r[i, QQQ[2]]+r[i, QQQ[3]], i = 1 .. n-1)])/sqrt(add(i, `in`(i, CovarianceMatrix(RB))));
LineChart(CumulativeSum([seq(r[i, QQQ[1]]+r[i, QQQ[2]]+r[i, QQQ[3]], i = 1 .. n-1)]), style = line, color = black, thickness = 2, labels = ["time", "Equity Curve Portfolio"], font = [Times, Roman, 14]);
 

########################################################

Ok, let me tell you what I am planing to do....

I was thinking I could have one slider that controlled all 5 element in the first correlation matrix. so for example if I increase 

the slider to 0.7 then all elemement in the correlation matrix would become 0.7 (except for the diagonal=1).

I would then apply Cholesky Decomposition to simulate 5*1000 random drawings which have a cross correlation defined by our

correlation matrix. I would then simulate 5 unit roots based upon the 5 columns of cross correlated random

drawings..... The plan is to have the 5 unot root plots updated when I change the cross correlation (one slider) 

To give you an example of what I have so far. (The code It is not working plus it only plot the first unit root but it is a start)

Not also that the slider should control the value of cor

restart;
with(LinearAlgebra);
with(Statistics);
C := Matrix(5, 5, symbol = corr, shape = symmetric);

myplot := proc (cor) local n, CD, A, B;

corr[1, 1] := 1; corr[2, 2] := 1; corr[3, 3] := 1; corr[4, 4] := 1; corr[5, 5] := 1;
corr[1, 2] := cor; corr[1, 3] := cor; corr[1, 4] := cor; corr[1, 5] := cor;
corr[2, 3] := cor; corr[2, 4] := cor; corr[2, 5] := cor; corr[3, 4] := cor;
corr[3, 5] := cor; corr[4, 5] := cor;
C;
IsDefinite(C);

with(LinearAlgebra);
with(ArrayTools);
with(Statistics);
randomize();
n := 1000;
CD := Matrix(LUDecomposition(evalf(C), 'method' = 'Cholesky'));
A := Matrix(Alias(Sample(RandomVariable(Normal(0, 3)), ColumnDimension(C)*n), [1 .. n, 1 .. ColumnDimension(C)]));
B := Matrix(Transpose(Multiply(CD, Transpose(A))));
CorrelationMatrix(B);

a1 := .2;
b := 1;
for i from 2 to n do
x[1] := 100;
x[i] := a1+b*x[i-1]+B[i, 1]
end do;

xx := [seq(x[i], i = 1 .. n)];
tt := [seq(i, i = 1 .. n)];
plot(tt, xx, color = black, thickness = 3, labels = ["time", "Stock 1"], font = [Times, Roman, 14])

end proc
 

Yes that worked great but

if I now want to simulate 1000 random drawings for each of the variables in the correlation matrix by using the code

with(LinearAlgebra);
with(ArrayTools);
with(Statistics);
randomize();

n := 10^3;
C := Matrix(5, 5, symbol = corr, shape = symmetric);
CD := Matrix(LUDecomposition(evalf(C), 'method' = 'Cholesky'), datatype = float[8]);
A := Matrix(Alias(Sample(RandomVariable(Normal(0, 3)), ColumnDimension(C)*n), [1 .. n, 1 .. ColumnDimension(C)]));
B := Matrix(Transpose(Multiply(CD, Transpose(A))));
CorrelationMatrix(B)

where should I put that code so that simulated random drawings are updated when I move the slider ??

ok, thanx Robert I will try that :-)

no this far excedes my ability. The only information I can extract is

gumbel = rstable(n, 1/theta) * (cos(pi/(2 * theta)))^theta

and somethink about a matrix and something about U and Y which I have no clue what that is.

Even if I spend the next 5 years trying to solve it ....I would probably not succed

I found I spreadsheet though that seams to work.... but it still uses a "black box" approach so it is quite usless if you want to really

understand what is going on....

Download 8342_Gumbel+Copula.zip
View file details

sorry to disappoint you but I am not sure I am skilled enough to be able to do that.

I found the section of code you are talking due to the good roadmap  but I can not make sense of it.

could you please explain in simple words what is going on?  :-) 

to tell you the truth I finding it a nightmare to install the QRMlib library.

It seams to be imposible to get it to work.....yet alone finding the copula function you are talking about.

I think I will just buy Matlab instead then I dont have to struggle to get this R stuff to work

How would this R copula look in Maple??

ok, yes but still if you call a R command like

"dcopula.gumbel(u, theta, logvalue=FALSE)"

what exactly is happening? hard to say?! (That I get a Gumbel copula is kind of given but how?)

R might be an option but I hate when I run in to something I dont understand.

I want to understand how it is done instead of simply using a pre programed black box function.

but maybe R or Matlabs black box commands is my only option if no one knows how it can be done ...hummmm

ok.

I think it is better to use copulas (even though I am not completly sure what it is exactly)

are you any good at that?  Do you know how to simulate correlated variables by using copulas?!

I know it can be done in Matlab...

what I mean is that if you run the below code then you will get a three variables that are correlated as we intialy have specified

(cor12=-0.6   cor13=-0.6  cor23=0.6)

restart;
with(LinearAlgebra);
with(ArrayTools);
with(Statistics);
n := 10^4;
cor12 := -.6; cor13 := -.6; cor23 := .6; cor21 := cor12; cor31 := cor13; cor32 := cor23;
C := Matrix([[1, cor12, cor13], [cor21, 1, cor23], [cor31, cor32, 1]]);
CD := Matrix(LUDecomposition(evalf(C), 'method' = 'Cholesky'), datatype = float[8]);
A := Matrix(Alias(Sample(RandomVariable(Normal(0, 1)), ColumnDimension(C)*n), [1 .. n, 1 .. ColumnDimension(C)]));
B := Multiply(CD, Transpose(A));
CorrelationMatrix(Transpose(B))

but if you change the correlation cor23 to -0.6 then we will get an error message saying

"Warning, Matrix is not positive-definite"

and your simulated correlation between 2 and 3 cor23 becomes horribly bad -0.37 

Thanx Alec !   :-)

Does it work for negative correlation as well (since we are using Cholesky decomp) (eigenvalue +)  ???

Parameter values for the Fractional Brownian Motion:

0 < H < 0.5        -  serial dependency

0.5 < H < 1        + serial dependency

H = 0.5               serial independent (pure random walk)

Yes I am definitely going to create some nice application around this one. (Maple Application Center) (with all credits due)

Your FBM implementation is to good to be ignored...