% CODE FOR TWO-SECTOR GROWTH MODEL

clear;

% hor = total number of periods
% T = periods before transition

hor=60;
T=20;


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

% Vectors, all containing zeros, ensuring they are equally long

Ge=zeros(1,hor);

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

YMx=zeros(1,hor);
YSx=zeros(1,hor);
Yx=zeros(1,hor);

Kx=zeros(1,hor);
Nx=zeros(1,hor);
Lx=zeros(1,hor);

wx=zeros(1,hor);
rx=zeros(1,hor);
Vx=zeros(1,hor);
zKx=zeros(1,hor);
zNx=zeros(1,hor);
nx=zeros(1,hor);
c1x=zeros(1,hor);
gKx=zeros(1,hor);

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

% "Free" parameters
mu=0.6;
phi=0.1;
beta=0.5;
gam=1.05;
L=1;

% Targeted parameters
nstar=gam^(1/(1-mu-phi));
wstar=1;
R=1;

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

% Initial conditions

KNaught = 0.0061;
NNaught = 0.0096;
AMNaught = 0.4320;
ASNaught = 0.4320;


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

% FUNCTIONS

% Factor shares

Vf=@(AM,AS,K,L) min(1,(((phi/(1-mu))^(1-mu))*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

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*(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*(K^(1-mu))*(N^(mu-1));

c1f=@(zK,zN,AM,K,N,L) (1-beta)*(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*(K^(phi-1))*(N^mu)*(L^(1-mu-phi));
 
% Function for updating K to next period

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

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

% Initializing state variables

AM=AMNaught;
AS=ASNaught; 
N=NNaught;
K=KNaught;

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

% Starting the loop

for t=1:1:hor;
    
    
    Ge(1,t)=t;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Step 1: compute variables for round t, save as scalars
    
    % Factor shares
    if Vf(AM,AS,K,L)<1;
    V=Vf(AM,AS,K,L);
    zN=zNf(AM,AS,K,L);
    zK=zKf(AM,AS,K,L);
    w=wMf(zK,zN,AM,K,N,L);
    % Nothing should change is calc w from Solow
    % w=wSf(zK,zN,AS,K,N);
    else
    V=1;
    zK=0;
    zN=0;
    w=wMf(zK,zN,AM,K,N,L);
    end
       
    % Output, consumption, capital growth, and interest rates
    YM=YMf(zK,zN,AM,K,N,L);
    YS=YSf(zK,zN,AS,K,N);
    c1=c1f(zK,zN,AM,K,N,L);
    gK=Kupdate(zK,zN,AM,K,N,L)/K;
    r=rMf(zK,zN,AM,K,N,L);
    n=nstar;

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Step 2: putting variables into vectors
    
    AMx(1,t)=AM;
    ASx(1,t)=AS;
    YMx(1,t)=YM;
    YSx(1,t)=YS;
    Yx(1,t)=YM+YS;
    zKx(1,t)=zK;
    zNx(1,t)=zN;
    Vx(1,t)=V;
    nx(1,t)=n;
    Nx(1,t)=N;
    Kx(1,t)=K;
    c1x(1,t)=c1;
    gKx(1,t)=gK;
    wx(1,t)=w;    
    rx(1,t)=r; 
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Step 3: updating K, N, AM, AS to next period
    
    K=Kupdate(zK,zN,AM,K,N,L);
    N=N*nx(1,t); 
    AM=AMNaught*gam^t;
    AS=ASNaught*gam^t;
    

end


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


% Figure
% Left as exercise