%
%    solow.m 
%    Economics 303:
%    Solow Model Graphs
%
%    George J Hall
%    Brandeis University
%    September 2006
%
clear
close all
%
%  set initial values for K, N and Y
%
K_1959 = 411.7;   % starting point for the capital stock
                  % net stock of nonresidential fixed private capital (in billions 1959 dollars)
                  % BEA Survey of Current Business, August 1994
                  % Fixed Reproducible Tangible Wealth in the United States
                    
N_1959 = 115.3;   % starting point for the labor force  
                  % 16+ population from Economic Report of the President 

Y_1959 = 506.6;   % starting point for output
                  % GDP (in billions of 1959 dollars)
                  % Economic Report of the President   
%
%  set values for parameters of the model
%
alpha = 1/3;   % capital share of income
delta = 0.08;  % depreciation rate of capital
s = 0.10;      % saving rate 
n = 0.015;     % population growth rate 
gamma = 0.0;   % growth rate in TFP
T = 100;       % number of periods to simulate 
%
%  back out initial level of A
%
A_1959 = Y_1959/( K_1959^(alpha)*N_1959^(1-alpha));
A_1959
%
%  initialize time series and 
%  compute initial capital-to-labor ratio
%
k = zeros(T,1);
N = zeros(T,1);
A = zeros(T,1);
k(1) = K_1959/N_1959;
N(1) = N_1959;
A(1) = A_1959; 
%
%  simulate model
%
for time = 1:T-1
    A(time+1,1) = (1+gamma)*A(time,1);
    N(time+1,1) = (1+n)*N(time,1);
    k(time+1,1) = (1-delta)*k(time,1)/(1+n) + s*A(time,1)*k(time,1)^(alpha)/(1+n);
end;
%
%  compute output per worker, total capital, total output,
%  wages and rental rate of capital
%
y = A.*k.^(alpha)*1000;
K = k.*N;
Y = A.*(K.^(alpha)).*(N.^(1-alpha));
rental_rate = alpha*A.*(K.^(alpha-1)).*(N.^(1-alpha));
wages = (1-alpha)*A.*(K.^(alpha)).*(N.^(-alpha));
%
%  plot a bunch of stuff
%
figure(1)
plot((1959:(1959+T-1)),k)
xlabel('year');
ylabel('capital-to-labor');
title('Simulated time path of the Capital-to-Labor Ratio');
print litk.ps

disp('pausing: hit the enter key');
pause

figure(2)
plot((1959:(1959+T-1)),y)
xlabel('year');
ylabel('output (1959 dollars)');
title('Simulated time path of real GDP per worker');
print lity.ps

disp('pausing: hit the enter key');
pause

figure(3)
plot((1959:(1959+T-1)),wages)
xlabel('year');
ylabel('wages (1959 dollars)');
title('Simulated time path of real wages');
print wages.ps

disp('pausing: hit the enter key');
pause

figure(4)
plot((1959:(1959+T-1)),rental_rate)
xlabel('year');
ylabel('rental rate of capital (1959 dollars)');
title('Simulated time path of real rental rate');
print rental_rate.ps

disp('pausing: hit the enter key');
pause

figure(5)
plot((1959:(1959+T-1)),Y)
xlabel('year');
ylabel('output (Billions of 1959 dollars)');
title('Simulated time path of real GDP');
print bigY.ps

disp('pausing: hit the enter key');
pause

figure(6)
plot((1959:(1959+T-1)),log(Y))
xlabel('year');
ylabel('log output');
title('Simulated time path of log real GDP');
print logbigY.ps

disp('pausing: hit the enter key');
pause

%
%  load in real GDP 
%
load GDPC1.DAT
T1 = length(GDPC1);

year = GDPC1(49:4:T1,1);
mth  = GDPC1(49:4:T1,2);
day  = GDPC1(49:4:T1,3);
RGDP = GDPC1(49:4:T1,4)*Y_1959/GDPC1(49,4); % Scale GDP to 1959 dollars. 
T2 = length(RGDP);

figure(7)
plot((1959:(1959+T2-1)),log(Y(1:T2)),year,log(RGDP),'r')
xlabel('year');
ylabel('log output');
title('Simulated and actual time path of log real GDP');
legend('simulated data','Scaled U.S. GDP',0); 
print logbigY.ps