%%%%%%%%%%%%%%%%%%%%  Solow Model (1956)  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
% Steady State.                                                           %
%                                                                         %
% [1] Relations of Solow Model (1956):                                    %
%   dk/dt = s*f(k(t))-theta*k(t)                                          %
%   netVest = s*f(k)-theta*k                                              %
%   k(0) = k_0 > 0                                                        %
%                                                                         %
% [2] C.E.S FUNCTION:                                                     %
% (C.E.S = Constant Elasticity of Substitution)                           %
%   Y = [a*K^(-gamma)+(1-a)*L^(-gamma)]^(-µ)                              %
%   with 0<a<1, gamma>-1 and different from zero.                         %
%                                                                         %
% The per capita formulation of the CES function                          %
% with constant returns to scale:                                         %
%             y = [a*k^(-gamma)+(1-a)]^(-1/gamma)                         %
%                                                                         %
% [3] COBB-DOUGLAS FUNCTION:                                              %
%   Y = A*K^(alpha)*L^(beta)                                              %
%   with alpha>0, beta<1.                                                 %
%   If beta = 1 - alpha, constant returns to scale. Then the (per capita) %
%   cobb-Douglas function is given by:                                    %
%                      y = A*k^(alpha)                                    % 
%                                                                         %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                                                                         %
%                                                                         %
% see also "ces.m", "cobbdouglas.m", "netVest.m", "TnetVest2.m", "Vest.m" %
%                                                                         %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


global A alpha gamma theta s;

 A=10; 
 alpha=0.32; 
 theta=0.09; % Depreciation rate 
 s=0.12; % Savings rate
 
% In the case gamma=0, we obtain a Cobb-Douglas function. 
% If gamma is different from 0, we obtain a CES function, and in such a
% case we suggest for example gamma=0.4.

disp(' ')
disp('We suggest gamma=0 to obtain a Cobb-Douglas function')
disp('... or gamma=0.4 for a CES function')
disp(' ')
gamma=input('gamma = ');
disp(' ')

% Find k* for which "netVest.m" is zero
kstar=fzero(@netVest,[10 80],0.0001)

% Savings, Investment.
k1=linspace(0,80,80);
Inv=Vest(k1);
dgn=theta*k1;
plot(k1,Inv,'b',k1,dgn,'r');
hold on;
Dstar=theta*kstar;
plot(kstar,Dstar,'r.','markersize',20);
xlabel(['Capital: k, where k* =', num2str(kstar)]);
ylabel('Invest.');
grid;
hold on;
lineStar=linspace(0,Dstar,30);
Qaxis=length(lineStar);
lineQ=ones(1,Qaxis).*kstar;
plot(lineQ,lineStar,'k-.');
hold off;

if gamma == 0; 
  title('Curves of "sf(k)" and "\thetak" where f(k) is a Cobb-Douglas function');
else
   title('Curves of "sf(k)" and "\thetak", where f(k) is a C.E.S function');
end;

disp(' ')
disp(' ')
disp('Press Enter To Continue ... ')
disp('.. to Stop the Simulation, Type "Ctrl-c" in the Command Window ..')
disp(' ')
disp(' ')

pause

% Steady state
% By using kstar=fzero(@netVest,[10 80],0.0001)

k=linspace(0,kstar,80);
NInv=netVest(k);
Nstar=netVest(kstar);

plot(k,NInv); 
xlabel(['Capital (k) with k* = ', num2str(kstar)]);
grid;
hold on;
plot(kstar,Nstar,'r.','markersize',20);
hold off;

if gamma == 0; 
  title('Net Investment: "sf(k)-\thetak", where f(k) is a Cobb-Douglas function');
else
   title('Net Investment: "sf(k)-\thetak", where f(k) is a C.E.S function');
end;

disp(' ')
disp('Press Enter To Continue ... ')
disp('.. to Stop the Simulation, Type "Ctrl-c" in the Command Window ..')
disp(' ')
disp(' ')

pause

 Ini_Stat=1;
 Tini=0;
 Tf=100;

[t,y]=ode45(@TnetVest2,[Tini Tf],Ini_Stat);

plot(t,y);
xlabel('Time');
ylabel('Capital');
grid;
hold on;
edy=y(end,:);
plot(Tf,edy,'k.','markersize',20);
hold off;

if gamma == 0; 
title('Time Evolution of net Invest. sf(k[t])-\thetak[t], when f(k) is a Cobb-Douglas function');
else
title('Time Evolution of net Invest. sf(k[t])-\thetak[t], when f(k) is a C.E.S function');  
end;


disp(' ')
disp('Press Enter To Continue ... ')
disp('.. to Stop the Simulation, Type "Ctrl-c" in the Command Window ..')
disp(' ')
disp(' ')

pause

 Ini_StatA =5;
 Ini_StatB =60;
 Tini=0;
 Tf=100; 
[tA,yA]=ode45(@TnetVest2,[Tini Tf],Ini_StatA);
[tB,yB]=ode45(@TnetVest2,[Tini Tf],Ini_StatB);
plot(tA,yA);
hold on;
plot(0,Ini_StatA,'b.','markersize',20);
hold on;
plot(0,Ini_StatB,'b.','markersize',20);
hold on;
plot(tB,yB);
xlabel('Time');
ylabel('Capital');
title('Time Evolution for 2 different initial values');
grid;
hold on;
edyA=yA(end,:);
edyB=yB(end,:);
plot(Tf,edyA,'k.','markersize',20);
hold on;
plot(Tf,edyB,'k.','markersize',20);
hold on;
lhA=length(tA);
lineA=ones(1,lhA).*kstar;
plot(tA,lineA,'k-.');
hold on;
lhB=length(tB);
lineB=ones(1,lhB).*kstar;
plot(tB,lineB,'k-.');
hold off


% Last revised on 10 March 2010.

%"Complex and Chaotic Nonlinear Dynamics. 
% Advances in Economics and Finance, 
% Mathematics and Statistics" 
% T.Vialar, Springer 2009.
% Copyright(c).