// \begin{verbatim}
// exec robust_diversification.sce

////////////////////////////////////////////
//  Land-use and diversification
////////////////////////////////////////////

function [y]=f(x,u,w)
// dynamics
  y=x*(1+u'*retu(w)); 
endfunction

function r=retu(w)
// return of land-use
  r=[r_0+sig*w;r_0-sig*w]; 
endfunction


xset("window",0);xbasc();
xtitle("Wealth",'time t','wealth x(t)');
// diversified
xset("window",11);xbasc();
xtitle("Wealth",'time t','wealth x(t)');
// specialized in use 1
xset("window",21);xbasc();
xtitle("Wealth",'time t','wealth x(t)');
// specialized in use 2

xset("window",2);xbasc();xtitle("Allocation",'t','u(t)');
xset("window",3);xbasc();xtitle("Uncertainty",'t','w(t)');

w_min=-1; w_max=1;
// uncertainty margins 

r_0=0.1;sig=0.5
// mean and variance return

Horizon=10; 
// time horizon

x_min=50; x_max=100;
// wealth bounds // x_min safety constraint 
r_star=r_0;
// Viable diversified return

x_prec=x_min/((1+r_star)^(Horizon-1));
// Viable initial state 

u_viab=[0.5;0.5]*ones(1,Horizon-1);
// Diversified viable strategy

u_1=[1;0]*ones(1,Horizon-1);
u_2=[0;1]*ones(1,Horizon-1);
// Specialized strategies

N_simu=10;
// Number of simulations

for i=1:N_simu
  // Simulations
  x_viab=x_prec+rand(1,Horizon)*(x_max-x_min);
  x_1=x_viab;
  x_2=x_viab;
  //  precautionnary initial conditions 
  // x(0)>= x_min/(1+r_star)^(Horizon-1)	
  w_viab=w_min+(w_max-w_min)*rand(1,Horizon-1);
  // Random climate along time
  for (t=1:1:Horizon-1)
	// Viable Trajectory x(.) u(.)
    x_viab(:,t+1)=f(x_viab(:,t),u_viab(:,t),w_viab(:,t));
    // Diversified wealth
    x_1(:,t+1)=f(x_1(:,t),u_1(:,t),w_viab(:,t));
    x_2(:,t+1)=f(x_2(:,t),u_2(:,t),w_viab(:,t));
	// Specialized wealth
  end
  //	
	rect1=[0,0,Horizon-1,200];
	rect2=[0,0,Horizon-2,1];
	rect3=[0,w_min,Horizon-2,w_max];
	xx=[0:1:Horizon-1];xu=[0:1:Horizon-2];
	//
	xset("window",0);
	plot2d(xx,[x_viab' x_min+zeros(1,Horizon)' ],rect=rect1);
	plot2d(xx,[x_viab' x_min+zeros(1,Horizon)' ],style=[1,-4])
	legends(['Diversified strategy';'Minimal wealth threshold'],...
			[1,-4],'ul')
	//
	xset("window",11);
	plot2d(xx,[x_1' x_min+zeros(1,Horizon)' ],rect=rect1) 
	plot2d(xx,[x_1' x_min+zeros(1,Horizon)' ],style=[1,-4])
	legends(['Strategy specialized in use 1';...
			 'Minimal wealth threshold'],[1,-4],'ul')
	//
	xset("window",21);
	plot2d(xx,[x_2' x_min+zeros(1,Horizon)' ],rect=rect1) 
	plot2d(xx,[x_2' x_min+zeros(1,Horizon)' ],style=[1,-4])
	legends(['Strategy specialized in use 2';...
			 'Minimal wealth threshold'],[1,-4],'ul')
	//
  xset("window",2);plot2d(xu,u_viab',rect=rect2);
  xset("window",3);plot2d(xu,w_viab',rect=rect3);
  
end
//\end{verbatim}