clear;

% hor = number of period per run
% cas = number of runs

hor=100;
cas=1000;


% Here I just define vectors containing zeros, ensuring they are equally long

Time=zeros(1,hor);

cx=zeros(cas,hor);
lx=zeros(cas,hor);
qx=zeros(cas,hor);
px=zeros(cas,hor);
rx=zeros(cas,hor);
Sx=zeros(cas,hor);

Diffx=zeros(cas,hor);

Countcx=zeros(cas,hor);
Countlx=zeros(cas,hor);

%parameters

z=0.2;

phi=0.5;

omega=1;

delta=0;
% Not used in paper; to allow competence to depreciate

% expressions

rhat=z/(phi*(1-z));
mu=(2*phi)*(1-(1-rhat)^.5)-(z/(1-z));

lthreshold=1-mu-z/(1-z);
% Level of l below which r=1, q=0, if c=1

Capomega=(2/rhat)*(1-(1-rhat)^.5)-1;

%Initial conditions

cNaught=.25;
lNaught=1.025*lthreshold

lNaught=0.7248;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Functions

r_thresh=z/(phi*(1-z));

rmax=@(c,l) ( z*(1-c)+(1-z)*(1-l) ) / (2*phi*(1-z));

cf=@(c,r,q) c*(1-delta)*(1-r)+r*q;

lf=@(l,r,q) l*(1-r)+r-phi*(r^2)*q;

Sf=@(c,l) omega*(z*c+(1-z)*l);

Psi=@(l) ((1-z)/z)*(1-mu-l);


% Starting the loop

for n=1:cas;

cx(n,1)=cNaught;
lx(n,1)=lNaught;


for v=1:hor-1;
    
     Time(v+1)=v;
     
     c=cx(n,v);
     l=lx(n,v);
     
     if c>.99;
     Countcx(n,v)=1;
     end
     
     if l>.99;
     Countlx(n,v)=1;
     end
     
     % Selecting some run where the transition happens in the middle
     if v<0.8*hor & v>0.6*hor & Countcx(n,v)==1;
         trans=n;
     end
     
     % Stage 3: determine if the incumbent survives
     Sx(n,v)=Sf(c,l);
     if rand>Sx(n,v);
         l=lNaught;
     end
     
    % Stage 4: update c, l; ruler chooses q, r
      if   c>Psi(l);
         q=1;
         r=rmax(c,l);
     else
         q=0;
         r=1;
     end
     
     qx(n,v)=q;
     rx(n,v)=r;
     px(n,v)=1-phi*r*q;
     
     Diffx(n,v)=c-Psi(l);
     
     cx(n,v+1)=cf(c,r,q);
     lx(n,v+1)=lf(l,r,q);
     
end

% Taking care of period "hor" 

qx(n,hor)=q;
rx(n,hor)=r;
px(n,hor)=1-phi*r*q;
Sx(n,hor)=Sf(cx(n,hor),lx(n,hor));

if cx(n,hor)>.99;
Countcx(n,hor)=1;
end 
if lx(n,hor)>.99;
Countlx(n,hor)=1;
end


end

     
SimFigure=subplot(2,2,1);

plot(Time,cx(trans,:),'b-*',Time,lx(trans,:),'r-+');
Leg1=legend('$c$','$l$');
legend('boxoff');
PeriodA=xlabel('Period');
ylim([-.1 1.1]);
PanA=ylabel('Fraction competent ($c$), loyal ($l$) in admin.');

subplot(2,2,2);
plot(Time,rx(trans,:),'b-*');
Leg2=legend('$r$');
legend('boxoff');
PeriodB=xlabel('Period');
ylim([-.1 1.1]);
PanB=ylabel('Fraction replaced ($r$)');

subplot(2,2,3);
plot(Time,Sx(trans,:),'b-*');
Leg3=legend('$S$');
legend('boxoff');
PeriodC=xlabel('Period');
ylim([-.1 1.1]);
PanC=ylabel('Survival probability ($S$)');

subplot(2,2,4);
plot(Time,px(trans,:),'b-*',Time,qx(trans,:),'r-+');
Leg4=legend('$p$','$q$');
legend('boxoff');
PeriodD=xlabel('Period');
ylim([-.1 1.1]);
PanD=ylabel('Fraction competent ($q$), loyal ($p$) recruits');


set(Leg1,'Interpreter','latex','FontSize',11,'FontName','times');
set(PanA,'Interpreter','latex','FontSize',11,'FontName','times');
set(Leg2,'Interpreter','latex','FontSize',11,'FontName','times');
set(PanB,'Interpreter','latex','FontSize',11,'FontName','times');
set(Leg3,'Interpreter','latex','FontSize',11,'FontName','times');
set(PanC,'Interpreter','latex','FontSize',11,'FontName','times');
set(Leg4,'Interpreter','latex','FontSize',11,'FontName','times');
set(PanD,'Interpreter','latex','FontSize',11,'FontName','times');

set(PeriodA,'Interpreter','latex','FontSize',11,'FontName','times');
set(PeriodB,'Interpreter','latex','FontSize',11,'FontName','times');
set(PeriodC,'Interpreter','latex','FontSize',11,'FontName','times');
set(PeriodD,'Interpreter','latex','FontSize',11,'FontName','times');

saveas(SimFigure,'C:\Work\Malthus2Stalin\MatLabFigs\SimFigure','fig');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ryiruy=lx(:,2);


% Monte Carlo stuff

meanc=mean(cx)';
meanl=mean(lx)';
meanq=mean(qx)';
meanp=mean(px)';
meanS=mean(Sx)';
meanr=mean(rx)';
meanCountc=mean(Countcx)';
meanCountl=mean(Countlx)';


MonteCarloFigure=subplot(2,2,1);
% 
plot(Time,meanc,'b-*',Time,meanl,'r-+');
LegN1=legend('Mean $c$','Mean $l$');
legend('boxoff');
PeriodNA=xlabel('Period');
ylim([-.1 1.1]);
PanNA=ylabel('Fraction competent ($c$), loyal ($l$) in admin.');

subplot(2,2,2);
plot(Time,meanCountc,'b-*',Time,meanCountl,'r-+');
LegN2=legend('Share $c>0.99$','Share $l>0.99$');
legend('boxoff');
PeriodNB=xlabel('Period');
ylim([-.1 1.1]);
PanNB=ylabel('Share of runs with $c>0.99$, $l>0.99$');

subplot(2,2,3);
plot(Time,meanS,'b-*',Time,meanr,'r-+');
LegN3=legend('Mean $S$','Mean $r$');
legend('boxoff');
PeriodNC=xlabel('Period');
ylim([-.1 1.1]);
PanNC=ylabel('Survival probability ($S$), fraction replaced ($r$)');

subplot(2,2,4);
plot(Time,meanp,'b-*',Time,meanq,'r-+');
LegN4=legend('Mean $p$','Mean $q$');
legend('boxoff');
PeriodND=xlabel('Period');
ylim([-.1 1.1]);
PanND=ylabel('Fraction competent ($q$), loyal ($p$) recruits');

set(LegN1,'Interpreter','latex','FontSize',11,'FontName','times');
set(LegN2,'Interpreter','latex','FontSize',11,'FontName','times');
set(LegN3,'Interpreter','latex','FontSize',11,'FontName','times');
set(LegN4,'Interpreter','latex','FontSize',11,'FontName','times');

set(PanNA,'Interpreter','latex','FontSize',11,'FontName','times');
set(PanNB,'Interpreter','latex','FontSize',11,'FontName','times');
set(PanNC,'Interpreter','latex','FontSize',11,'FontName','times');
set(PanND,'Interpreter','latex','FontSize',11,'FontName','times');

set(PeriodNA,'Interpreter','latex','FontSize',11,'FontName','times');
set(PeriodNB,'Interpreter','latex','FontSize',11,'FontName','times');
set(PeriodNC,'Interpreter','latex','FontSize',11,'FontName','times');
set(PeriodND,'Interpreter','latex','FontSize',11,'FontName','times');

saveas(MonteCarloFigure,'C:\Work\Malthus2Stalin\MatLabFigs\MonteCarloFigure','fig');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Sensitivity analysis
% Change parameter values above to get the time paths in the plot below

SensitivityAnalysisFigure=subplot(2,2,1);
 
plot(Time,meanc,'b-*');
%LegS1=legend('Baseline','$z=0.15$','$\phi=0.95$','$\omega=0.95$');
PeriodSA=xlabel('Period');
ylim([-.1 1.1]);
PanSA=ylabel('Mean level of $c$');

subplot(2,2,2);
plot(Time,meanl,'b-*');
%LegS2=legend('Mean $c$','Mean $S$','Mean $r$');
PeriodSB=xlabel('Period');
ylim([-.1 1.1]);
PanSB=ylabel('Mean level of $l$');

subplot(2,2,3);
plot(Time,meanS,'b-*');
%LegS3=legend('Mean $c$','Mean $S$','Mean $r$');
PeriodSC=xlabel('Period');
ylim([-.1 1.1]);
PanSC=ylabel('Mean level of $S$');

subplot(2,2,4);
plot(Time,meanr,'b-*');
%LegS4=legend('Mean $c$','Mean $S$','Mean $r$');
PeriodSD=xlabel('Period');
ylim([-.1 1.1]);
PanSD=ylabel('Mean level of $r$');

%set(LegS1,'Interpreter','latex','FontSize',11,'FontName','times');

set(PanSA,'Interpreter','latex','FontSize',11,'FontName','times');
set(PanSB,'Interpreter','latex','FontSize',11,'FontName','times');
set(PanSC,'Interpreter','latex','FontSize',11,'FontName','times');
set(PanSD,'Interpreter','latex','FontSize',11,'FontName','times');

set(PeriodSA,'Interpreter','latex','FontSize',11,'FontName','times');
set(PeriodSB,'Interpreter','latex','FontSize',11,'FontName','times');
set(PeriodSC,'Interpreter','latex','FontSize',11,'FontName','times');
set(PeriodSD,'Interpreter','latex','FontSize',11,'FontName','times');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





















