/* 課題4　解答例 */

/* ［数値処理］→［自動的に数値で出力］*/
if numer#false then numer:false else numer:true;
/*浮動小数点の桁数=7に指定 */
fpprec : 7;
/* 浮動小数点表示桁数を7に指定*/
fpprintprec:7;

/* (I) */
/* 期待収益率を定義 */
mu:[0.1, 0.05, 0.15]; 
/* 分散共分散行列を定義 */
V : matrix([0.04,-0.015,0.01],[-0.015,0.01,-0.02],[0.01,-0.02,0.09]);
Rf:0.04;

m: floor((lmax(mu)-lmin(mu))/0.001) 
/* 束縛条件の mu_p を muの[最小値，最大値]で0.001 刻みで設定 
floor(x)はx以下の最大整数を返す*/ ;
mu_p: makelist(lmin(mu)+i*0.001, i, 0, m);
if mu_p[m] < lmax(mu) then mu_p: endcons(lmax(mu), mu_p);
/* muの(最小値-最大値)が0.01で割り切れなかった場合の処理 */
m:length(mu_p);

var_p:makelist(0,m) /* 最小分散リスト初期化 */;
weight: zeromatrix(m,3) /* 最小分散投資比率リスト初期化*/;

/* 不等式制約付き最小化パッケージfmim_coblya読み込み */
load(fmin_cobyla);
for i:1 thru m do (
result: fmin_cobyla(
[w1,w2,w3].V.transpose([w1,w2,w3]), 
[w1,w2,w3], 
[0.3,0.3,0.4], 
constraints=[
[1,1, 1].[w1,w2,w3]=1, 
mu.[w1,w2,w3]=mu_p[i]
]), 
weight[i]: result[1], var_p[i]: result[2]
);

sigma_p: sqrt(var_p); weight;

sigma_mu:makelist([sigma_p[i], mu_p[i]], i, 1, length(mu_p));
slope:lmax((mu_p-Rf)/sigma_p);
f(x):=Rf+slope*x;
load(draw);
draw2d(
points(sigma_mu), 
explicit(f(x), x, 0, lmax(sigma_p)), 
yrange=[0.04, 0.16], 
xlabel="sigma", ylabel="mu", 
title="Efficient frontier and minimum variance boundary (3 risky asset portfolio)"
);

/* (II) */
load(fmin_cobyla);
for i:1 thru m do (
result: fmin_cobyla(
[w1,w2,w3].V.transpose([w1,w2,w3]), 
[w1,w2,w3], 
[0.3,0.3,0.4], 
constraints=[
[1,1, 1].[w1,w2,w3]=1, 
mu.[w1,w2,w3]=mu_p[i], 
w1>=0, w2>=0, w3>=0
]), 
weight[i]: result[1], var_p[i]: result[2]
);

sigma_p: sqrt(var_p); weight;

sigma_mu:makelist([sigma_p[i], mu_p[i]], i, 1, length(mu_p));
slope:lmax((mu_p-Rf)/sigma_p);
f(x):=Rf+slope*x;
load(draw);
draw2d(
points(sigma_mu), 
explicit(f(x), x, 0, lmax(sigma_p)), 
yrange=[0.04, 0.16], 
xlabel="sigma", ylabel="mu", 
title="Efficient frontier and minimum variance boundary (3 risky asset portfolio; no shortsales)"
);