clc             % clear screen
clear all       % clear all variables
 
%
% parameter values
%
% natural survival rates 
b=exp(-[0,0.2,0.2,0.2,0.2,0.2,0.2,0.2]);         
% relative prices at age
p=ones(size(b));
p(1)=0;
% average price per kg of fish
paverage=1;
% weights at age
w=[0.0,0.166,0.409,0.783,1.189,1.704,2.578,3.826]; % in kg
% age-specific fecundities = weight x fraction mature
m=[0.13,0.36,0.83,0.94,0.96,0.96,0.96,0.98];
% age-specific catchabilities
phi=[0.0,0.6,0.8,1.0,1.0,1.0,1.0,1.0];
% parameters of smooth hockey-stick stock-recruitment function
a=[420.038,1.642/420.038];
%R(SSB) = Rmax (1-exp(-alpha/Rmax * SSB)
% cost parameter
catchability=6.57e-6;       % 
c=.554e-3/catchability;     % 0.554 tons/day, so cost parameter is in 1000 tons
% initial stock
x0=[173.198,173.198,191.178,132.910,67.820,21.887,11.748,4.000]; % 1000 individuals 2010
ssbmsy=997; % BMSY_cons, Froese/Proelss (2010)
msy=255; % MSY_cons, Froese/Proelss (2010)
% time horizon from 2010 onwards
T=31;
%
 
%
% Simulations
%
 
%
% 15%-Option
%
x=zeros(8,T); % stock numbers at age
x(:,1)=x0;    
harvest=56.8;
F=zeros(T,1); % fishing mortality
for t=2:T
 if harvest < (1-exp(-0.3))*w*(phi'.*x(:,t-1))
     F(t-1)=harvest/(w*(phi'.*x(:,t-1)));
     harvest=1.15*harvest;
 else
     F(t-1)=1-exp(-0.3);
 end
 x(1,t) = a(1)*(1-exp(-a(2)*(w.*m)*x(:,t-1)));  % 1-year-old: smooth hockey stick SSR 
 for age=2:7
  x(age,t)=b(age-1)*(1-F(t-1)*phi(age-1))*x(age-1,t-1);              
  end
  x(8,t)=b(7)*(1-F(t-1)*phi(7))*x(7,t-1)+b(8)*(1-F(t-1)*phi(8))*x(8,t-1); % 8+-year-old    
end
if harvest < (1-exp(-0.3))*w*(phi'.*x(:,T))
    F(T)=harvest/(w*(phi'.*x(:,T)));    
else
    F(T)=1-exp(-0.3);     
end
ssb_15percent=((w.*m)*x)';    
harvest_15percent=F.*(w.*phi*x)';
profit_15percent=paverage*(F.*(p.*w.*phi*x)'+c*log(ones(T,1)-F));  
 
%
% F03-Option
%
x=zeros(8,T);    
x(:,1)=x0;
F=(1-exp(-0.3))*ones(T,1);
 
F(1)=56.8/(w*(phi'.*x0'));
for t=2:T       
    x(1,t) = a(1)*(1-exp(-a(2)*(w.*m)*x(:,t-1)));   % 1-year-old: smooth hockey stick SSR 
    for age=2:7
        x(age,t) = b(age-1)*(1-F(t-1)*phi(age-1))*x(age-1,t-1);              
    end    
    x(8,t) = b(7)*(1-F(t-1)*phi(7))*x(7,t-1)+b(8)*(1-F(t-1)*phi(8))*x(8,t-1); % 8+-year-old    
end    
ssb_F03=((w.*m)*x)';    
harvest_F03=F.*(w.*phi*x)';
profit_F03=paverage*(F.*(p.*w.*phi*x)'+c*log(ones(T,1)-F));  
 
%
% F06-Option
%
x=zeros(8,T);    
x(:,1)=x0;
F=(1-exp(-0.6))*ones(T,1);
 
F(1)=56.8/(w*(phi'.*x0'));
for t=2:T       
    x(1,t) = a(1)*(1-exp(-a(2)*(w.*m)*x(:,t-1)));   % 1-year-old: smooth hockey stick SSR 
    for age=2:7
        x(age,t) = b(age-1)*(1-F(t-1)*phi(age-1))*x(age-1,t-1);              
    end    
    x(8,t) = b(7)*(1-F(t-1)*phi(7))*x(7,t-1)+b(8)*(1-F(t-1)*phi(8))*x(8,t-1); % 8+-year-old    
end    
ssb_F06=((w.*m)*x)';    
harvest_F06=F.*(w.*phi*x)';
profit_F06=paverage*(F.*(p.*w.*phi*x)'+c*log(ones(T,1)-F));  
 
 
%
% HCR-Option
%
x=zeros(8,T); % stock numbers at age
x(:,1)=x0;    
hmin=56.8;
F=zeros(T,1); % fishing mortality
for t=2:T
   ssb=(w.*m)*x(:,t-1);
   if ssb<0.5*ssbmsy
        harvest=hmin;
   elseif ssb>ssbmsy
        harvest=0.9*msy;
   else
       harvest=hmin+(0.9*msy-hmin)*(ssb-0.5*ssbmsy)/(0.5*ssbmsy);
   end   
 F(t-1)=harvest/(w*(phi'.*x(:,t-1)));
 
 x(1,t)=a(1)*(1-exp(-a(2)*ssb));  % 1-year-old: smooth hockey stick SSR 
 for age=2:7
  x(age,t)=b(age-1)*(1-F(t-1)*phi(age-1))*x(age-1,t-1);              
  end
  x(8,t)=b(7)*(1-F(t-1)*phi(7))*x(7,t-1)+b(8)*(1-F(t-1)*phi(8))*x(8,t-1); % 8+-year-old    
end
   ssb=(w.*m)*x(:,T);
   if ssb<0.5*ssbmsy
        harvest=hmin;
   elseif ssb>ssbmsy
        harvest=0.9*msy;
   else
       harvest=hmin+(0.9*msy-hmin)*(ssb-0.5*ssbmsy)/(0.5*ssbmsy);
   end   
 F(T)=harvest/(w*(phi'.*x(:,T)));
 
ssb_hcr=((w.*m)*x)';
harvest_hcr=F.*(w.*phi*x)';
harvest_hcr(T)=harvest;
profit_hcr=paverage*(F.*(p.*w.*phi*x)'+c*log(ones(T,1)-F));  
 
%
% write and plot output
%
output=[(2010:2010+T-1)',...    
    ssb_F03(1:T),harvest_F03(1:T),profit_F03(1:T),...
    ssb_15percent(1:T),harvest_15percent(1:T),profit_15percent(1:T),...
    ssb_hcr(1:T),harvest_hcr(1:T),profit_hcr(1:T),...    
    ssb_F06(1:T),harvest_F06(1:T),profit_F06(1:T)];
 
save('Output_2011_02_11.mat','output');
xlswrite('Output_2011_02_11.xls',output);
 
plot((2010:2010+T-1),harvest_15percent,'-x',...
    (2010:2010+T-1),harvest_F03,'-x',...
    (2010:2010+T-1),harvest_hcr,'-x',...
    (2010:2010+T-1),harvest_F06,'-x')
pause
plot((2010:2010+T-1),ssb_15percent,'-x',...
    (2010:2010+T-1),ssb_F03,'-x',...
    (2010:2010+T-1),ssb_hcr,'-x',...
    (2010:2010+T-1),ssb_F06,'-x')
pause
plot((2010:2010+T-1),profit_15percent,'-x',...
    (2010:2010+T-1),profit_F03,'-x',...
    (2010:2010+T-1),profit_hcr,'-x',...
    (2010:2010+T-1),profit_F06,'-x')