% NOTE: SOME FIGURES ARE SAVED AS .FIG FILES BELOW, SO YOU HAVE TO NAME
% FOLDERS ACCORDINGLY ON THE MACHINE WHERE YOU RUN THIS CODE


clear;

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

hor=14;
cas=500;


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

Ge=zeros(1,hor);

War=ones(cas,hor);

AMx=zeros(cas,hor);
ASx=zeros(cas,hor);

YM1x=zeros(cas,hor);
YS1x=zeros(cas,hor);
Y1x=zeros(cas,hor);

YM2x=zeros(cas,hor);
YS2x=zeros(cas,hor);
Y2x=zeros(cas,hor);

K1x=zeros(cas,hor);
N1x=zeros(cas,hor);
L1x=zeros(cas,hor);

K2x=zeros(cas,hor);
N2x=zeros(cas,hor);
L2x=zeros(cas,hor);

w1x=zeros(cas,hor);
r1x=zeros(cas,hor);
V1x=zeros(cas,hor);
zK1x=zeros(cas,hor);
zN1x=zeros(cas,hor);
n1x=zeros(cas,hor);
c11x=zeros(cas,hor);
gK1x=zeros(cas,hor);

w2x=zeros(cas,hor);
r2x=zeros(cas,hor);
V2x=zeros(cas,hor);
zK2x=zeros(cas,hor);
zN2x=zeros(cas,hor);
n2x=zeros(cas,hor);
c12x=zeros(cas,hor);
gK2x=zeros(cas,hor);

xunknownx=zeros(cas,hor);
xknownx=zeros(cas,hor);

dK1x=zeros(cas,hor);
dK2x=zeros(cas,hor);
dN1x=zeros(cas,hor);
dN2x=zeros(cas,hor);


%parameters

mu=0.6;
phi=0.1;
gamM=1.032;
gamS=1.518;
% To eliminate transition, set gamS=1;
beta=0.5;
tao=0.01;

%War related parameters

%Mweight same as theta in paper; also misnomer
Mweight=1;

% zeta=1 means dK=0.3 and dN=0.06 (half of "hat" values)
zeta=1;
dKhat=.6; 
dNhat=.12;

%Omega = risk aversion parameter in Gov'ts' utility
omega=0;

% T is last period when zK and zN are zero in no-war simulation
T=5;

%n*=fertility on Malthusian BGP.
nstar=gamM^(1/(1-mu-phi));

% w*=wage on Malthusian BGP; normalized to one
wstar=1;


%Initial conditions
AMNaught=1;
LNaught=.5;

%Setting initial N and K so that economy is on a Malthusian BGP: n*=gK* and
%initial wage equals w*
R=mu/(beta*(1-phi));
NNaught=LNaught*((mu*(1-tao)*AMNaught/wstar)*(wstar/(R*nstar))^phi)^(1/(1-mu-phi));
KNaught=(wstar/(R*nstar))*NNaught;

% Initial AS is set so that Solow sector is not active at first; becomes
% active in period T+1 in no-war simulation

ASNaught=(gamS^(-T))*((phi/(1-mu))^(1-mu))*(1/((1-tao)*mu))*((R*nstar)^(1-mu))*(wstar)^mu;

%FUNCTIONS

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

%Factor shares

Expr=(phi/(1-mu))^(1-mu);

Vf=@(AM,AS,K,L) min(1,(Expr*AM/AS)^(1/(1-phi-mu))*(L/K));

zKf=@(AM,AS,K,L) 1-Vf(AM,AS,K,L);

zNfprel=@(zK,V) zK/(zK + V*(1-mu)/phi);

zNf=@(AM,AS,K,L) zNfprel(zKf(AM,AS,K,L),Vf(AM,AS,K,L));

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

% Incomes, wages, and consumption

% Note: wM and wS are equal, but when zK and zN are (close to) one (V close to zero) we cannot
% compute wages from wM, thus two functions; this matters when time spans are long (e.g. 100 periods)

YMf=@(zK,zN,AM,K,N,L) ((1-zK)^phi)*((1-zN)^mu)*AM*(K^phi)*(N^mu)*(L^(1-mu-phi));

wMf=@(zK,zN,AM,K,N,L) mu*((1-zK)^phi)*((1-zN)^(mu-1))*AM*(1-tao)*(K^phi)*(N^(mu-1))*(L^(1-mu-phi));

YSf=@(zK,zN,AS,K,N) (zK^(1-mu))*(zN^mu)*AS*(K^(1-mu))*(N^mu);

wSf=@(zK,zN,AS,K,N) mu*(zK^(1-mu))*(zN^(mu-1))*AS*(1-tao)*(K^(1-mu))*(N^(mu-1));

c1f=@(zK,zN,AM,K,N,L) (1-beta)*(1-tao)*(1/N)*YMf(zK,zN,AM,K,N,L)*(1-phi-mu+mu/(1-zN));

% Interest rate, computed from the Malthus sector

rMf=@(zK,zN,AM,K,N,L) phi*((1-zK)^(phi-1))*((1-zN)^mu)*AM*(1-tao)*(K^(phi-1))*(N^mu)*(L^(1-mu-phi));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% Computing war gains

% Land conquest function: Conq
% Same as H/2 in paper

Conq=@(K1,K2,N1,N2) ((K1/(K1+K2))^(1-mu))*((N1/(N1+N2)))^mu;

% Possible alterations of Conq for sens analysis

% To create perm peace:
% Conq=@(K1,K2,N1,N2) 0;

% To make conquests indep. of K, N; set Conq=1/2 (or H=1 in paper), i.e.,
% as if symmetric strengths

% Conq=@(K1,K2,N1,N2) 1/2;

% Conq=@(K1,K2,N1,N2) ((K1^(1-mu))*(N1^mu))/((K1^(1-mu))*(N1^mu)+(K2^(1-mu))*(N2^mu));

% Death and depreciation functions

dKf=@(K1,K2,N1,N2) (dKhat/2)*(zeta+(1-zeta)*2*Conq(K1,K2,N1,N2));

dNf=@(K1,K2,N1,N2) (dNhat/2)*(zeta+(1-zeta)*2*Conq(K1,K2,N1,N2));

% Factor allocations with war

zKwarf=@(AM,AS,K,L,ta,dK) zKf(AM,AS,K*(1-dK),L+ta*(1-L));

zNwarf=@(AM,AS,K,L,ta,dK) zNf(AM,AS,K*(1-dK),L+ta*(1-L));

%  Malthusian output in war

YMwarf=@(AM,AS,K,N,L,ta,dK,dN) YMf(zKwarf(AM,AS,K,L,ta,dK),zNwarf(AM,AS,K,L,ta,dK),AM,K*(1-dK),N*(1-dN),L+ta*(1-L));

%  Solowian output in war

YSwarf=@(AM,AS,K,N,L,ta,dK,dN) YSf(zKwarf(AM,AS,K,L,ta,dK),zNwarf(AM,AS,K,L,ta,dK),AS,K*(1-dK),N*(1-dN));

% Utility function for government

U=@(G) (1/(1-omega))*(G^(1-omega));


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

% Updating K to next period

Kupdate=@(zK,zN,AM,K,N,L) beta*(1-tao)*YMf(zK,zN,AM,K,N,L)*(1-phi-mu+mu/(1-zN));

%Function for fertility behavior

eta=1.5;
nuw=7;
nbar=1.5;

b=(eta*nstar-nbar)/(eta-1);
a=(nbar-1)/(nuw-eta);
d=(nuw*nbar-eta)/(nuw-eta);

nwf1=@(w) (nstar-b)*(w/wstar)+b;
nwf2=@(w) d-a*(w/wstar);

% Starting the loop

for m=1:cas;
    
    %Initializing state variables before starting new run
    
    opp=0;
    
    AM=AMNaught;
    AS=ASNaught;
    
    N1=NNaught;
    K1=KNaught;
    L1=LNaught;
    
    N2=NNaught;
    K2=KNaught;
    L2=LNaught;
  

for v=1:hor;
    
    
    %The only state variables that differ btw the 2 countries are K, L, and N. 
      
      L1x(m,v)=L1;
      L2x(m,v)=L2;
      
      % WAR PHASE
     
      opp=rand;
      % xknown is distr equally on 0, 0.1, ..., 0.9
      xknown = -.1+ceil((12-1).*rand(1,1))./10;
      % To take out known shock; must set higher than mean of known x
      % (=1/2), to make sure war breaks out 1st period
      % xknown=0.6;
      xknownx(m,v)=xknown;
     
     % Death and depr rates in war
     % Note 1st argument = K2 for country 1, since it fights against 2's K and N.
     % Same way 3rd argument is N2 for country 1 etc.
     
        dK1=dKf(K2,K1,N2,N1); 
        dN1=dNf(K2,K1,N2,N1); 
        dK2=dKf(K1,K2,N1,N2); 
        dN2=dNf(K1,K2,N1,N2); 
        
        dK1x(m,v)=dK1;
        dK2x(m,v)=dK2;
        dN1x(m,v)=dN1;
        dN2x(m,v)=dN2;
    
% COMPUTING POTENTIAL WAR GAINS IN sow STATES OF THE WORLD

sow=11;

% Vectors with outcomes in all sow states of the world

YM1_x=zeros(sow,1);
YM2_x=zeros(sow,1);
YS1_x=zeros(sow,1);
YS2_x=zeros(sow,1);

YMwar1x=zeros(sow,1);
YMwar2x=zeros(sow,1);
YSwar1x=zeros(sow,1);
YSwar2x=zeros(sow,1);

Ugain1x=zeros(sow,1);
Ugain2x=zeros(sow,1);

 for c=1:sow;
     
     % The unknown shock takes values 0, 0.1,...,1 with equal prob
     
     xunknownC=((c-1)/(sow-1));
     
     % ta1 and ta2 same as g in paper
     
     % Country 1
     
     ta1=min(.99,xknown*xunknownC*2*Conq(K1,K2,N1,N2));
     
     YMwar1x(c,1)=YMwarf(AM,AS,K1,N1,L1,ta1,dK1,dN1);
     YM1_x(c,1)=YMf(zKf(AM,AS,K1,L1),zNf(AM,AS,K1,L1),AM,K1,N1,L1);
     
     if Vf(AM,AS,K1,L1)<1;
        YSwar1x(c,1)=YSwarf(AM,AS,K1,N1,L1,ta1,dK1,dN1);
        YS1_x(c,1)=YSf(zKf(AM,AS,K1,L1),zNf(AM,AS,K1,L1),AS,K1,N1);
     else
        YSwar1x(c,1)=0;
        YS1_x(c,1)=0;
     end
     
     Ugain1x(c,1)=U(YMwar1x(c,1)+Mweight*YSwar1x(c,1))-U(YM1_x(c,1)+Mweight*YS1_x(c,1));
     
     % Country 2
     
     ta2=min(.99,xknown*xunknownC*2*Conq(K2,K1,N2,N1));
     
     YMwar2x(c,1)=YMwarf(AM,AS,K2,N2,L2,ta2,dK2,dN2);
     YM2_x(c,1)=YMf(zKf(AM,AS,K2,L2),zNf(AM,AS,K2,L2),AM,K2,N2,L2);
     
     if Vf(AM,AS,K2,L2)<1;
        YSwar2x(c,1)=YSwarf(AM,AS,K2,N2,L2,ta2,dK2,dN2);
        YS2_x(c,1)=YSf(zKf(AM,AS,K2,L2),zNf(AM,AS,K2,L2),AS,K2,N2);  
     else
        YSwar2x(c,1)=0;
        YS2_x(c,1)=0;
     end
     
     Ugain2x(c,1)=U(YMwar2x(c,1)+Mweight*YSwar2x(c,1))-U(YM2_x(c,1)+Mweight*YS2_x(c,1));
     
 end
 

 
 % Expected war gains
 
Gain1=mean(Ugain1x(:,1));
Gain2=mean(Ugain2x(:,1));
     
     % Determining actual conquests by letting the unknown shock be
     % realized
     % xunknown is distr equally on 0, 0.1, ..., 0.9
     xunknown = -.1+ceil((12-1).*rand(1,1))./10;
     xunknownx(m,v)=xunknown;
     
     ta1=min(.99,xknown*xunknown*2*Conq(K1,K2,N1,N2));
     ta2=min(.99,xknown*xunknown*2*Conq(K2,K1,N2,N1));
     
     % Making the war decision
      
     if opp>1/2 & Gain1>0;
         L1x(m,v)=L1x(m,v)+ta1*L2x(m,v);
         L2x(m,v)=L2x(m,v)*(1-ta1);
         K1=(1-dK1)*K1;
         K2=(1-dK2)*K2;
         N1=(1-dN1)*N1;
         N2=(1-dN2)*N2;
     elseif opp<1/2 & Gain2>0;
         L2x(m,v)=L2x(m,v)+ta2*L1x(m,v);
         L1x(m,v)=L1x(m,v)*(1-ta2);
         K1=(1-dK1)*K1;
         K2=(1-dK2)*K2;
         N1=(1-dN1)*N1;
         N2=(1-dN2)*N2;
     else
         War(m,v)=0;
     end
     
    
    L1=L1x(m,v);
    L2=L2x(m,v);
    
    
    Ge(1,v)=v;
    
    % PRODUCATION PHASE
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Factor shares and wages
    
    %Remembering that zK and zN are constrained to zero if V = (constrained to) 1.
    
    % Country 1
    
    if Vf(AM,AS,K1,L1)<1;
    V1=Vf(AM,AS,K1,L1);
    zN1=zNf(AM,AS,K1,L1);
    zK1=zKf(AM,AS,K1,L1);
    w1=wMf(zK1,zN1,AM,K1,N1,L1);
    % Nothing should change is calc w1 from Solow
    % w1=wSf(zK1,zN1,AS,K1,N1);
    else
    V1=1;
    zK1=0;
    zN1=0;
    w1=wMf(zK1,zN1,AM,K1,N1,L1);
    end
    
    % Country 2
    
    if Vf(AM,AS,K2,L2)<1;
    V2=Vf(AM,AS,K2,L2);
    zN2=zNf(AM,AS,K2,L2);
    zK2=zKf(AM,AS,K2,L2);
    w2=wMf(zK2,zN2,AM,K2,N2,L2);
    % Nothing should change is calc w2 from Solow
    % w2=wSf(zK2,zN2,AS,K2,N2);
    else
    V2=1;
    zK2=0;
    zN2=0;
    w2=wMf(zK2,zN2,AM,K2,N2,L2);
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Output, consumption, capital growth, and interest rates
    
    % Country 1

    YM1=YMf(zK1,zN1,AM,K1,N1,L1);
    YS1=YSf(zK1,zN1,AS,K1,N1);
    c11=c1f(zK1,zN1,AM,K1,N1,L1);
    gK1=Kupdate(zK1,zN1,AM,K1,N1,L1)/K1;
    r1=rMf(zK1,zN1,AM,K1,N1,L1);
    
    % Country 2
    
    YM2=YMf(zK2,zN2,AM,K2,N2,L2);
    YS2=YSf(zK2,zN2,AS,K2,N2);
    c12=c1f(zK1,zN2,AM,K2,N2,L2);
    gK2=Kupdate(zK2,zN2,AM,K2,N2,L2)/K2;
    r2=rMf(zK2,zN2,AM,K2,N2,L2);
    
    % Fertility
    
    % Country 1
    
    if w1<eta*wstar;
        n1=nwf1(w1);
    elseif w1>eta*wstar & w1<nuw*wstar;
        n1=nwf2(w1);
    else
        n1=1;
    end
    
    % Country 2
    
    if w2<eta*wstar;
        n2=nwf1(w2);
    elseif w2>eta*wstar & w2<nuw*wstar;
        n2=nwf2(w2);
    else
        n2=1;
    end
    
    AMx(m,v)=AM;
    ASx(m,v)=AS;
    
    % Putting variables into vectors
    
    % Country 1
    
    YM1x(m,v)=YM1;
    YS1x(m,v)=YS1;
    Y1x(m,v)=YM1+YS1;
    zK1x(m,v)=zK1;
    zN1x(m,v)=zN1;
    V1x(m,v)=V1;
    n1x(m,v)=n1;
    N1x(m,v)=N1;
    K1x(m,v)=K1;
    c11x(m,v)=c11;
    gK1x(m,v)=gK1;
    w1x(m,v)=w1;    
    r1x(m,v)=r1; 
    
    % Country 2
    
    YM2x(m,v)=YM2;
    YS2x(m,v)=YS2;
    Y2x(m,v)=YM2+YS2;
    zK2x(m,v)=zK2;
    zN2x(m,v)=zN2;
    V2x(m,v)=V2;
    n2x(m,v)=n2;
    N2x(m,v)=N2;
    K2x(m,v)=K2;
    c12x(m,v)=c12;
    gK2x(m,v)=gK2;
    w2x(m,v)=w2;   
    r2x(m,v)=r2;
    
    %Updating K and N to next period
    
    % Country 1
    
    K1=Kupdate(zK1,zN1,AM,K1,N1,L1);
    N1=N1*n1;
    
    % Country 2
    
    K2=Kupdate(zK2,zN2,AM,K2,N2,L2);
    N2=N2*n2;
    
    %Updating AM and AL to next period
    
    AM=AMNaught*gamM^v;
    AS=ASNaught*gamS^v;
    

end

end


%Fertility net of war deaths

n1xw=n1x.*(1-War.*dN1);
n2xw=n2x.*(1-War.*dN2);

% Change definition of periods

% In a no-war simulation, period T+1 is the last period in which zK is 
% constrained to 0
% Explainer: Ge starts at 1 and zK>0 exactly T+1 periods after period 1.
% That is, in period T+2. T+1 is thus the last period in which zK=0
% We now make make 0 the last period in in which zK is constrained to 0 
% in the no-war simulation

Ge=Ge-(T+1);

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

%Figure A1: one single run (here the very first)
% Fig 2 in paper
 
FigureA1=subplot(2,2,1);

plot(Ge,n1x(1,:),'b-*',Ge,n1xw(1,:),'b:*',Ge,n2x(1,:),'r-',Ge,n2xw(1,:),'r:', 'LineWidth',2)
A1subh1=legend('Fertility, country 1', 'Pop. growth, country 1','Fertility, country 2','Pop. growth, country 2');
legend('boxoff');
xlim([-5 7]);
A1subr1=title('Panel A: fertility ($n$) and population growth');

subplot(2,2,2);
plot(Ge,L1x(1,:),'b-',Ge,L2x(1,:),'r:','LineWidth',2);
A1subh2=legend('Country 1', 'Country 2');
legend('boxoff');
xlim([-5 7]);
A1subr2=title('Panel B: territorial holdings');

subplot(2,2,3);
semilogy(Ge,w1x(1,:),'b-',Ge,w2x(1,:),'r:',Ge,eta*wstar*ones(1,hor),'k',Ge,nuw*wstar*ones(1,hor),'k','LineWidth',2);
A1subh3=legend('Country 1','Country 2');
legend('boxoff');
xlim([-5 7]);
A1subr3=title('Panel C: wages ($w$)');

subplot(2,2,4);
plot(Ge,zK1x(1,:),'b-*',Ge,zN1x(1,:),'b:*',Ge,zK2x(1,:),'r-',Ge,zN2x(1,:),'r:');
A1subh4=legend('$z_{K}$, country 1','$z_{N}$, country 1','$z_{K}$, country 2','$z_{N}$, country 2',2);
legend('boxoff');
xlim([-5 7]);
A1subr4=title('Panel D: shares of capital and labor in Solow sector');

set(A1subh1,'Interpreter','latex');
set(A1subh2,'Interpreter','latex');
set(A1subh3,'Interpreter','latex');
set(A1subh4,'Interpreter','latex');

set(A1subr1,'Interpreter','latex');
set(A1subr2,'Interpreter','latex');
set(A1subr3,'Interpreter','latex');
set(A1subr4,'Interpreter','latex');


saveas(FigureA1,'C:\Work\War\RandRSummer2009\MatlabFigs\FigureA1','fig');

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

% AGGREGATING (taking the average across runs)

% War frequency

Warfreq=mean(War);

% Selecting which country wins/loses in each run

for m=1:cas;
if L1x(m,hor)>L2x(m,hor);
    
    nwinner(m,:)=n1x(m,:);
    nloser(m,:)=n2x(m,:);
    nwinnerw(m,:)=n1xw(m,:);
    nloserw(m,:)=n2xw(m,:);
    
    Lwinner(m,:)=L1x(m,:);
    
    Nwinner(m,:)=N1x(m,:);
    
    Wwinner(m,:)=w1x(m,:);
    Wloser(m,:)=w2x(m,:);
    
    ZKwinner(m,:)=zK1x(m,:);
    ZKloser(m,:)=zK2x(m,:);
    
    rwinner(m,:)=r1x(m,:);
    rloser(m,:)=r2x(m,:);
    
    Ywinner(m,:)=Y1x(m,:);
    Yloser(m,:)=Y2x(m,:);
    
    dKwinner(m,:)=dK1x(m,:);
    dKloser(m,:)=dK2x(m,:);
    dNwinner(m,:)=dN1x(m,:);
    dNloser(m,:)=dN2x(m,:);
    
else
    nwinner(m,:)=n2x(m,:);
    nloser(m,:)=n1x(m,:);
    nwinnerw(m,:)=n2xw(m,:);
    nloserw(m,:)=n1xw(m,:);
    
    Lwinner(m,:)=L2x(m,:);
    
    Nwinner(m,:)=N2x(m,:);
    
    Wwinner(m,:)=w2x(m,:);
    Wloser(m,:)=w1x(m,:);
    
    ZKwinner(m,:)=zK2x(m,:);
    ZKloser(m,:)=zK1x(m,:);
    
    rwinner(m,:)=r2x(m,:);
    rloser(m,:)=r1x(m,:);
    
    Ywinner(m,:)=Y2x(m,:);
    Yloser(m,:)=Y1x(m,:);
    
    dKwinner(m,:)=dK2x(m,:);
    dKloser(m,:)=dK1x(m,:);
    dNwinner(m,:)=dN2x(m,:);
    dNloser(m,:)=dN1x(m,:);
    
end
end

nwinnerag=mean(nwinner)';
nloserag=mean(nloser)';
nwinneragw=mean(nwinnerw)';
nloseragw=mean(nloserw)';
Lwinnerag=mean(Lwinner)';
Nwinnerag=mean(Nwinner)';
Lloserag=1-Lwinnerag;
ZKwinnerag=mean(ZKwinner)';
ZKloserag=mean(ZKloser)';
Wwinnerag=mean(Wwinner)';
Wloserag=mean(Wloser)';
rwinnerag=mean(rwinner)';
rloserag=mean(rloser)';
Ywinnerag=mean(Ywinner)';
Yloserag=mean(Yloser)';
dKwinnerag=mean(dKwinner)';
dKloserag=mean(dKloser)';
dNwinnerag=mean(dNwinner)';
dNloserag=mean(dNloser)';

% Done aggregating

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

% Figure A2: similar to A1 but average across cas runs
% Fig 3 in paper

FigureA2=subplot(2,2,1);
plot(Ge,nwinnerag,'b-*',Ge,nwinneragw,'b:*',Ge,nloserag,'r-',Ge,nloseragw,'r:','LineWidth',2)
xlim([-5 7]);
ylim([0.9 1.5]);
A2subh1=legend('Fertility, winner', 'Pop. growth, winner','Fertility, loser','Pop. growth, loser');
legend('boxoff');
A2subr1=title('Panel A: fertility ($n$) and population growth');

subplot(2,2,2);
plot(Ge,Lwinnerag,'b',Ge,Lloserag,'r','LineWidth',2);
xlim([-5 7]);
A2subh2=legend('Winner', 'Loser');
legend('boxoff');
A2subr2=title('Panel B: territorial holdings');

subplot(2,2,3);
semilogy(Ge,Wwinnerag,'b-',Ge,Wloserag,'r:',Ge,eta*ones(1,hor),'k',Ge,nuw*ones(1,hor),'k','LineWidth',2);
A2subh3=legend('Winner','Loser');
xlim([-5 7]);
legend('boxoff')
A2subr3=title('Panel C: wages ($w$)');

subplot(2,2,4);
plot(Ge,Warfreq,'k-*',Ge,ZKwinnerag,'b-',Ge,ZKloserag,'r:','LineWidth',2);
A2subh4=legend('Fraction at war','$z_{K}$, winner','$z_{K}$, loser',2);
xlim([-5 7]);
legend('boxoff');
A2subr4=title('Panel D: fraction at war and share of capital in Solow sector');

set(A2subh1,'Interpreter','latex');
set(A2subh2,'Interpreter','latex');
set(A2subh3,'Interpreter','latex');
set(A2subh4,'Interpreter','latex');

set(A2subr1,'Interpreter','latex');
set(A2subr2,'Interpreter','latex');
set(A2subr3,'Interpreter','latex');
set(A2subr4,'Interpreter','latex');

saveas(FigureA2,'C:\Work\War\RandRSummer2009\MatlabFigs\FigureA2','fig');

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

%Figure A3: land/population shares, and p.c. income gaps
% Fig 4 in paper

% Note: last period with zK=0 is period T+1

hor3=T+1+5;
hor2=T+1;
hor1=1;

N1Share1=N1x(:,hor1)./(N1x(:,hor1)+N2x(:,hor1));
N1Share2=N1x(:,hor2)./(N1x(:,hor2)+N2x(:,hor2));
N1Share3=N1x(:,hor3)./(N1x(:,hor3)+N2x(:,hor3));

PCShareIII=(Y1x(:,hor3)./N1x(:,hor3))./(Y2x(:,hor3)./N2x(:,hor3));
PCShareII=(Y1x(:,hor2)./N1x(:,hor2))./(Y2x(:,hor2)./N2x(:,hor2));
PCShareI=(Y1x(:,hor1)./N1x(:,hor1))./(Y2x(:,hor1)./N2x(:,hor1));

Shares=horzcat(N1Share3,N1Share2,N1Share1);

FigureA3=subplot(1,2,1);

hist(Shares,11);
h1=title('Histograms over population share of country 1','FontSize',12);
A3subh1=legend('5 periods after transition','Period of transition','First period','FontSize',12);
set(h1,'Interpreter','latex');
set(A3subh1,'Interpreter','latex');
legend('boxoff');

subplot(1,2,2);
plot(L1x(:,hor3),PCShareIII,'rx',L1x(:,hor2),PCShareII,'go',L1x(:,hor1),PCShareI,'b*','MarkerSize',10);
A4subt1=ylabel('Per-capita income gap','FontSize',12);
A4subt2=xlabel('Land share','FontSize',12);
A4subh1=title('Land share and relative per-capita income of country 1','FontSize',12);
A4subr1=legend('5 periods after transition','Period of transition','First period','FontSize',12);
legend('boxoff');

set(A4subh1,'Interpreter','latex');
set(A4subr1,'Interpreter','latex');
set(A4subt1,'Interpreter','latex');
set(A4subt2,'Interpreter','latex');

saveas(FigureA3,'C:\Work\War\RandRSummer2009\MatlabFigs\FigureA3','fig');






















