exec("utils.sci");// charge les primitives graphiques utiles

// On définit les caracteristiques des actifs
d=2;
mu=[0.05,0.15];
// Matrice de covariance: des 1 sur la diagonale, des rho ailleurs
rho=0.0;
covariance=rho*ones(d,d)+(1-rho)*eye(d,d);
sigma=[0.10,0.30];
Gamma = diag(sigma) * covariance * diag(sigma);

// Les caractéristiques des actifs de base
moyenne_actif=mu;
std_actif=sigma;

// Tracé des points représentant les actifs dans le plan 
//                 (ecart-type,moyenne)
rectangle=plot_actifs(moyenne_actif,std_actif);
// On recupére le "rectangle" dans lequel on va dessiner par la suite.
//-------------------------------------------------------
// On parcours toutes les valeurs possibles de x_1 et x_2
// C'est facile lorsque d=2 ...
N=100; 
x_1=[0:1/N:1];
x=[x_1;1-x_1];

moyenne_x=0;std_x=0;
for i=[1:N+1] do
  current_x = x(:,i);
  QUESTION: moyenne_x(i) = QUE VAUT LA MOYENNE DE CE PORTEFEUILLE
  QUESTION: std_x(i)= QUE VAUT L'ECART-TYPE DE CE PORTEFEUILLE
end;

plot2d(std_x, moyenne_x, style=1,rect=rectangle);
set_line(bleu,1);// on trace la courbe en bleu, épaisseur=1
//-------------------------------------------------------
// Calcul du point de variance minimum
[m,imin]=min(std_x);

plot2d(std_x(imin), moyenne_x(imin), rect=rectangle);// trace du point 
set_dot(vert,4);// en vert, taille=4
//-------------------------------------------------------
moyenne_x=0;std_x=0;
// On autorise l'emprunt de l'actif 1
//  -> la quantité d'actif 1 est négative

N=1000;

// QUESTION: générer N x_1 entre -5 et 0 puis les x_2 correspondants;  
x_1=[-5:A_COMPLETER:0];
x_2=A_COMPLETER;
x   = [x_1;x_2];

for i=[1:N+1] do
  current_x = x(:,i);
  moyenne_x(i+1)=mu * current_x;
  std_x(i+1)=sqrt(current_x'*Gamma*current_x);
end;

plot2d(std_x, moyenne_x,style=point, rect=rectangle);
set_line(vert,1);// trace les lignes en vert, taille=1
//-------------------------------------------------------
// Que se passe t'il lorsque l'on emprunte de l'actif 2 ?
moyenne_x=0;std_x=0;

// On autorise l'emprunt de l'actif 2
//  -> la quantité d'actif 2 est négative
N=1000;
// QUESTION: générer N x_2 entre -5 et 0 puis les x_1 correspondants;  
x_2 = [-5:A_COMPLETER:0];
x_1 = A_COMPLETER;
x   = [x_1;x_2];

for i=[1:N+1] do
  current_x = x(:,i);
  moyenne_x(i+1)=mu * current_x;
  std_x(i+1)=sqrt(current_x'*Gamma*current_x);
end;

plot2d(std_x, moyenne_x,style=point, rect=rectangle);
set_line(rouge,1);// trace les lignes en rouge, taille=1
//-------------------------------------------------------
// Matrice de covariance: des 1 sur la diagonale, des rho ailleurs
rho=0.5;// -0.5
covariance=rho*ones(d,d)+(1-rho)*eye(d,d);
Gamma = diag(sigma) * covariance * diag(sigma);
// etc ...
//-------------------------------------------------------
// On rajoute un actif sans risque de moyenne nulle
r0=0;
mu=[r0,mu];
// comme ce rendement est suppose deterministe, la matrice de  
// variance covariance se complete par une ligne et une colonne de 0
QUESTION: Gamma= QUE VAUT GAMMA DANS CE CAS

moyenne_actif=mu;
std_actif=sqrt(diag(Gamma));

// On materialise les 3 actifs de base
plot2d(std_actif, moyenne_actif, style=1, rect=rectangle);
set_dot(noir,4);

// On considère des portefeuilles *avec l'actif sans risque*
// mais *sans emprunt*. On les tire au hasard dans le simplexe
// de dimension 3
N=1000;
moyenne_x=0;std_x=0;
for i=[1:N]
  x=simplexe(d+1);// tirage au hasard dans le simplexe
  moyenne_x(i)=mu * x;
  std_x(i)=sqrt(x'*Gamma*x);
end;

plot2d(std_x, moyenne_x,style=-1, rect=rectangle);
set_dot(bleuclair,2);
//-------------------------------------------------------
// On genere des portefeuilles sans actifs sans risque
N=1000;
moyenne_y=0;std_y=0;
for i=[1:N]  
  QUESTION: y = LES PORTEFEUILLES SANS ACTIF SANS RISQUE

  // on calcule les moyennes et variances des portefeuilles y
  moyenne_y(i)=mu * y;
  std_y(i)=sqrt(y'*Gamma*y);
end;

// Le point P maximise la pente de la droite entre (sigma0=0, x_0=r0) 
// et les portefeuilles y (sans actif sans risque)
r0=mu(1);
sigma0=0;// sigma0 = Gamma(1,1) = 0
pente=(moyenne_y - r0) ./ (std_y - sigma0);// calcul des pentes
[lambda,imax]=max(pente);

// Tracé du point P
x_P=moyenne_y(imax);
sigma_P=std_y(imax);
plot2d(sigma_P,x_P, style=1, rect=rectangle);
set_dot(vert,4);

// Tracé du segment "Actif sans risque -> P"
plot2d([sigma0,sigma_P],[r0,x_P], style=1, rect=rectangle);
set_line(vert,2);

// Tracé de la droite "actif sans risque -> P" au dela de P
lambda=(x_P-r0)/(sigma_P-sigma0);//pente de la droite
sigma=2.0;// arbitraire mais "grand"
x=r0+lambda*(sigma-sigma0);
plot2d([sigma0,sigma],[r0,x], rect=rectangle);
set_line(vert,2);
//-------------------------------------------------------
// L'emprunt en actif sans risque est autorisé
N=1000;
moyenne_x=0;variance_x=0;
for i=[1:N]
  // On génère des portefeuilles dont la quantité 
  // d'actif sans risque est uniforme sur [-4,1]
  x_0=grand(1,1,'unf',-4,1);
  s=simplexe(d);// tirage uniforme dans le simplexe de dim |$d$|.
  // On veut génèrer un portefeuille avec x_0 actifs sans risque
  // et de valeur totale |$x_0 + \sum_{i=1,\ldots,d} x_i = 1$|.
  QUESTION: x = COMMENT DEFINIR CE PORTEFEUILLE EN UTILISANT s

  moyenne_x(i)=mu * x;
  variance_x(i)=sqrt(x'*Gamma*x);
end;

// Tracé des points tirés au hazard
plot2d(variance_x, moyenne_x,style=-2,rect=rectangle);
set_dot(mauve,2);
//-------------------------------------------------------
d=3;
N=1000;
moyenne_x=0;variance_x=0;
for i=[1:N]
  x=arbitraire(d);
  moyenne_x(i)=mu * x;
  variance_x(i)=sqrt(x'*Gamma*x);
end;

// Tracé des points tirés au hazard
plot2d(variance_x, moyenne_x,style=-2,rect=rectangle);
set_dot(mauve,2);
//-------------------------------------------------------
// Caracteristiques des actifs sans risque

d=3;rho=0.0;

min_esp=0.05;max_esp=0.15;
mu=[min_esp:(max_esp-min_esp)/(d-1):max_esp];

// On suppose que tous les actifs risqués ont une 
// corrélation constante égale à |$\rho$|.
// On doit forcement avoir |$\rho >= -(1/(d-1))$|, 
// sinon la matrice n'est pas une matrice de covariance (exercice!).
covariance= rho*ones(d,d)+(1-rho)*eye(d,d);

// On choisit un ecart type croissant en fonction de l'actif
min_sigma=0.1;max_sigma=0.3;
sigma=[min_sigma:(max_sigma-min_sigma)/(d-1):max_sigma];

// La matrice de variance covariance se calcule par :
Gamma = diag(sigma) * covariance * diag(sigma);

// Les caractéristiques des actifs de base
moyenne_actif=mu;
std_actif=sqrt(diag(Gamma));

// Tracé des actifs dans le plan (ecart-type,moyenne)
rectangle=plot_actifs(moyenne_actif,std_actif);
//-------------------------------------------------------
// On simule des valeurs possibles 
N=1000;
moyenne_x=0;std_x=0;
for i=[1:N]
  x=simplexe(d);
  moyenne_x(i)=mu * x;
  std_x(i)=sqrt(x'*Gamma*x);
end;

plot2d(std_x, moyenne_x, style=1,rect=rectangle);
set_dot(bleu,2);

// Le point de variance minimum
[m,imin]=min(std_x);
plot2d(std_x(imin), moyenne_x(imin), rect=rectangle);
set_dot(vert,4);
//-------------------------------------------------------
N=1000;
moyenne_x=0;std_x=0;
sigma_e=2;
for i=[1:N] do
  x=arbitraire(d,sigma_e);
  moyenne_x(i)=mu * x;
  std_x(i)=sqrt(x'*Gamma*x);
end;

plot2d(std_x, moyenne_x,style=point, rect=rectangle);
f=gcf();Dessin=f.children(1).children(1).children(1);
set_dot(rouge,2);
//-------------------------------------------------------
r0=0;
mu=[r0,mu];
// comme le rendement est deterministe, la matrice de 
// variance covariance se complete par
Gamma=[zeros(1,d+1);zeros(1,d)',Gamma];

// On materialise les 4 actifs de base
moyenne_actif=mu;
std_actif=sqrt(diag(Gamma));
plot2d(std_actif, moyenne_actif, style=1, rect=rectangle);
set_dot(noir,4);
//-------------------------------------------------------
N=1000;
std_x=0;moyenne_x=0;
for i=[1:N]
  x=simplexe(d+1);
  moyenne_x(i)=mu * x;
  std_x(i)=sqrt(x'*Gamma*x);
end;

plot2d(std_x, moyenne_x,style=-1, rect=rectangle);
set_dot(bleuclair,2);
//-------------------------------------------------------
N=1000;
moyenne_y=0;std_y=0;
for i=[1:N]
  y = [0;simplexe(d)];    // un portefeuille sans actif sans risque
  moyenne_y(i) = mu * y;  // on calcule ses caractéristiques
  std_y(i)=sqrt(y'*Gamma*y);
end;

r0=mu(1);
sigma0=sqrt(Gamma(1,1));// sigma0 = Gamma(1,1) = 0
pente=(moyenne_y - r0) ./ (std_y - sigma0);// calcul des pentes
[lambda,imax]=max(pente);
x_P=moyenne_y(imax);
sigma_P=std_y(imax);
//-------------------------------------------------------
// Tracé du point P
plot2d(sigma_P,x_P, style=1, rect=rectangle);
set_dot(red,2);

// Tracé du segment "Actif sans risque -> P"
plot2d([sigma0,sigma_P],[r0,x_P], style=1, rect=rectangle);
set_line(vert,2);

// Tracé de la droite "actif sans risque -> P" au dela de P
sigma=2.0;// arbitraire mais "grand"
lambda=(x_P-r0)/(sigma_P-sigma0);//pente de la droite
x=r0+lambda*(sigma-sigma0);
plot2d([sigma0,sigma],[r0,x],style=-1, rect=rectangle);
set_line(vert,2);
//-------------------------------------------------------
N=5000;
moyenne_x=0;variance_x=0;
for i=[1:N]
  x=arbitraire(d+1,sigma_e);
  moyenne_x(i)=mu * x;
  variance_x(i)=sqrt(x'*Gamma*x);
end;

// Tracé des points tirés au hazard
plot2d(variance_x, moyenne_x,style=-2,rect=rectangle);
set_dot(mauve,2);
//-------------------------------------------------------
// Choix des caracteristiques des actifs sans risque
d=30;rho=0.0;
min_esp=0.05;max_esp=0.15;
mu = [min_esp:(max_esp-min_esp)/(d-1):max_esp];
min_sigma=0.1;max_sigma=0.3;
sigma = [min_sigma:(max_sigma-min_sigma)/(d-1):max_sigma];
correlation = rho*ones(d,d)+(1-rho)*eye(d,d);
Gamma = diag(sigma) * correlation * diag(sigma);

// Définition de la fonction de cout et de sa dérivée
// en utilisant la contrainte |$\sum_{i=1}^d \lambda_i=1$|
function [f,g,ind]=cost(x,ind)
// On maximise 
//       f = (mu*lambda)^2 / lambda'*Gamma*lambda
// sous la contrainte |$\sum_{i=1}^d \lambda_i=1$|, |$x=\lambda(2:d)$|
  p=prod(size(x))+1;// dimension de |$\lambda = 1+dim(x)$|
  // x_1 est calculé en fonction de |$\lambda$| à partir de x(2:d)
  // en utilisant la contrainte |$\sum_{i=1}^d \lambda_i = 1$|
  lambda=[1-sum(x);x];
  ps=mu*lambda;
  var=lambda'*Gamma*lambda;
  // On cherche à minimiser (lambda.mu)^2 / lambda'.Gamma.lambda
  f=ps^2 / var;
  // le dérivée du coût en fonction de lambda
  k=(2*ps/var)*mu - (2*ps^2/var^2)*lambda'*Gamma;
  // Calcul de la dérivée par rapport a |$x$| 
  // en fonction de la derivee en |$\lambda$|.
  g=k(2:p)-k(1);
  f=-f;g=-g;// On maximise mais Scicoslab suppose que l'on minimise ...
endfunction

x0=ones(d-1,1)/d;
[f,xopt]=optim(cost,x0);
Xopt=[1-sum(xopt);xopt];
Fopt=sqrt(-f)

// Est on bien entre |$0$| et |$1$| ? C'est toujours le cas pour 
// Rho diagonale mais pas toujours dans le cas non diagonal
ok = and(0<=Xopt) & and(Xopt<=1)

// Lorsque rho=0, il y a une solution explicite (exercice)
// lambda_i = alpha * mu_i/sigma_i^2, renormalisé
// On verifie ...
sigma=sqrt(diag(Gamma))';
x=(mu ./ sigma^2);
x=x'/sum(x);// noramlisation
norm(x - Xopt) // lorsque la matrice Rho est diagonale 
               // ca devrait etre petit
//-------------------------------------------------------
function lambda=x2lambda(x,R,mu)
  // construction de lambda a partir de x
  // en tenant compte des contraintes |$\sum_i \lambda_i  =1$|, |$r^T \lambda = R$|
  p=prod(size(x))+2;// dimension de lambda = 2+dim x
  alpha0=1-sum(x);
  beta0=R-mu(1:p-2)*x;
  lambda_n_1 = (mu(d)*alpha0 - beta0) / (mu(d) - mu(d-1));
  lambda_n = (beta0 - mu(d-1)*alpha0) / (mu(d) - mu(d-1));
  lambda=[x;lambda_n_1;lambda_n];
endfunction

function [f,g,ind]=cost(x,ind)
// On minimise la variance 
//       f = lambda'*Gamma*lambda
// sous les contraintes |$\sum(\lambda)=1$|, |$r^T \lambda = R$|
  p=prod(size(x))+2;// dimension de lambda = 2+dim x
  lambda=x2lambda(x,R,mu);
  k=Gamma*lambda;// derive en fonction de lambda
  f=lambda'*k;// = lambda'*Gamma*lambda;
  for i=[1:p-2] do
    // derivee par rapport a x (et non lambda)
    g(i)=k(i) ...
      + (k(p-1) * (- mu(p) + mu(i)) + k(p) * (- mu(i) + mu(p-1))) ...
          / (mu(d) - mu(d-1));
  end    
endfunction

// Choix des caracteristiques des actifs sans risque
d=30;rho=0.0;
min_esp=0.05;max_esp=0.15;
mu=[min_esp:(max_esp-min_esp)/(d-1):max_esp];
min_sigma=0.1;max_sigma=0.3;
sigma=[min_sigma:(max_sigma-min_sigma)/(d-1):max_sigma];
correlation = rho*ones(d,d)+(1-rho)*eye(d,d);
Gamma = diag(sigma) * correlation * diag(sigma);

x0=ones(d-2,1)/d;
i=0;abscisse=0;ordonnee=0;
for R=[0.05:0.001:0.15] do
  [f,xopt]=optim(cost,x0);
  Xopt=x2lambda(xopt,R,mu);
  i=i+1;
  abscisse(i)=sqrt(f);
  ordonnee(i)=R;
end;
plot2d(abscisse,ordonnee);
//-------------------------------------------------------