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


//parameters
puis=0.5;
Horizon=10;
rho=0.95 ;
B0=10;

p=0.1; q=1-p;
p=0.9; q=1-p;

R#=1.3;Rb=1;

// isoelastic utility 
function [u]=util(x)
  u=x^puis
endfunction

// certainty equivalent
Rchap=(p*util(R#)+q*util(Rb))^(1/puis) ;

a=(rho * util(Rchap))^(1/(puis-1)) ;

// proportion of consumption  b(t)
b=[];
b(Horizon+1)=1;
for t=Horizon:-1:1
  b(t)=a*b(t+1)/(1+a*b(t+1)) ;
end

// optimal catches and resource
function [h]=PREL_OPT(t,B)
  h =b(t)*B
endfunction

function [x]=DYN_OPT(t,B,R)
  x=R*(B-PREL_OPT(t,B))
endfunction

xset("window",1) ; xbasc()
xset("window",2) ; xbasc()


N_simu=3;
for i=1:N_simu

// simulation of random productivity R(t)
z=rand(1,Horizon,'uniform'); 
for t=1:Horizon
  if (z(t) <= p) then 
	Rsimu(t)=R#;
  else Rsimu(t)=Rb;
  end
end

// computation of optimal trajectories 
Bopt(1)=B0;J=0;
for t=1: Horizon
  hopt(t)=PREL_OPT(t,Bopt(t));
  Bopt(t+1)=DYN_OPT(t,Bopt(t),Rsimu(t));
end

// graphics
xset("window",1) ; 
plot2d(1:Horizon+1,Bopt,style=i) ;
plot2d(1:Horizon+1,Bopt,style=-i) ;
xtitle("Optimal biomass","time t","B(t)")

xset("window",2) ; 
plot2d(1:Horizon,hopt,style=i) ;
plot2d(1:Horizon,hopt,style=-i) ;
xtitle("Optimal catch","time t","h(t)")

end


// \end{verbatim}