Maple Example

maple_ex.html

>	with(linalg);

>	M:=matrix(3,4,[60,50,53,28,110,120,127,59,30,30,30,15]);

>	gausselim(M);

>	gaussjord(M);

>	P:=Vector[row]([5/10,1/10,4/10]);

>	m:=[(6/25)/P[1],(3/50)/P[2],(1/5)/P[3]];

>	R1:=Vector[row]([M[1,1]/M[1,4],M[1,2]/M[1,4],M[1,3]/M[1,4]]);

>	R2:=Vector[row]([M[2,1]/M[2,4],M[2,2]/M[2,4],M[2,3]/M[2,4]]);

>	R3:=Vector[row]([M[3,1]/M[3,4],M[3,2]/M[3,4],M[3,3]/M[3,4]]);

>	ER1:=dotprod(R1,P);#E(R1)

>

evalf(%);

>	ER2:=dotprod(R2,P);#E(R2)

>

evalf(%);

>	dotprod(R3,P);

>	s1:=sqrt(P[1]*(R1[1]-ER1)^2+(P[2])*(R1[2]-ER1)^2+(P[3])*(R1[3]-ER1)^2);#stdev(R1)

>

evalf(%);

>	s2:=sqrt(P[1]*(R2[1]-ER2)^2+(P[2])*(R2[2]-ER2)^2+(P[3])*(R2[3]-ER2)^2);#stdev(R2)

>

evalf(%);

>	R1R2:=[R1[1]*R2[1],R1[2]*R2[2],R1[3]*R2[3]];

>	ER1R2:=dotprod(R1R2,P);#E(R1R2)

>	(ER1R2-ER1*ER2)/(s1*s2);#rho(R1,R2)

>

evalf(%);

>	mR1:=[m[1]*R1[1],m[2]*R1[2],m[3]*R1[3]];

>	dotprod(mR1,P);#check that E(mR1)=1

>	mR2:=[m[1]*R2[1],m[2]*R2[2],m[3]*R2[3]];

>	dotprod(mR2,P); #check that E(mR2)=1

>	Em:=dotprod(m,P); #E(m)=1/Rf

>	sm:=sqrt(P[1]*(m[1]-Em)^2+P[2]*(m[2]-Em)^2+P[3]*(m[3]-Em)^2);#stdev(m)

>

evalf(%);

>	(1-Em*ER1)/(sm*s1);#rho(m,R1)

>

evalf(%);

>	(1-Em*ER2)/(sm*s2);#rho(m,R2)

>

evalf(%);

>	yu:=1/Em+sm/Em*x;#upper part of MV frontier yu=Rf+sigma(m)/E(m)*x

>	yl:=1/Em-sm/Em*x;#lower part of MV frontier yl=Rf-sigma(m)/E(m)*x

>	Ex:=ER1*w+(1-w)*ER2;#E(Rp), Rp=wR1+(1-w)R2

>	sigma:=sqrt(w^2*s1^2+(1-w)^2*s2^2+2*w*(1-w)*(ER1R2-ER1*ER2));#stdev(Rp)

>	plot([[x,yu(x),x=0..0.2],[x,yl(x),x=0..0.2],[sigma(w),Ex(w),w=0..1]],color=[red,red,blue]);

>	msq:=[(12/25)^2, (3/5)^2, (1/2)^2];#calculate m^2

>	dotprod(msq,P);#E(m^2)

>

evalf(%);

>	moverEmsq:=625/157*[12/25, 3/5, 1/2];

>	dotprod(moverEmsq,P);#E(m/E(m^2))

>

evalf(%);

>	mxmoverEmsq:=[m[1]*moverEmsq[1],m[2]*moverEmsq[2],m[3]*moverEmsq[3]];

>	dotprod(mxmoverEmsq,P);#check E(m*m/E(m^2))=1

>	v1:=[1,0,-6/5];

>	u2:=[1,-4,0];

>	v1u2:=[1,0,0];

>	v1sq := [1, 0, 36/25];

>	v2 := u2-dotprod(v1u2, P)*v1/dotprod(v1sq, P);

>	v1v2 := [144/269, -4*0, -6/5*150/269];

>	dotprod(P, v1v2);

>	v2sq := [20736/72361, 16, 22500/72361];

>	Restar := dotprod(P, v1)*v1/dotprod(P, v1sq)+dotprod(P, v2)*v2/dotprod(P, v2sq);

>	plot3d(-1/4*x2-5/6*x3, x2 = (-10 .. 10), x3 = (-10 .. 10),axes=NORMAL);

>	[2, 2, 2] = moverEmsq+ga*Restar;

>	g := (2-300/157)/(7/157);

>	(2-375/157)/(-61/314);

>	(2-625/314)/(3/628);

>	dst:=diff(sigma,w);

>	dE:=diff(Ex,w);

>	dEoverdst := proc (w) options operator, arrow; 174*sqrt(4898091*w^2+1269296-4848088*w)/(9796182*w-4848088) end proc;

>	plot({dEoverdst(w),sqrt(3)/25},w=.5..1);

>	solve(dEoverdst(w)=sqrt(3)/25);

>

evalf(%);

>