% This program implements the calculations in class notes 11, of Theory of Income, Fall 2004, UofC, see file notes11theoryincomefall2004.tex
% Written by F Alvarez

clear all; close all;

% adds the result of the current simulation to the file pset7.txt
% next line deletes the content of prevoius simulations 
! del neolabor.txt 
diary('neolabor.txt'); 

% Parameters (all variables are in yearly units)
% **********

% Utility function: v(c,n) = ( c^(1-sigma) - 1 )/(1-sigma) - A n^(1+phi)/(1+phi) ;
sigma = 1  ;   % risk aversion
phi   = 0 ;  % elasticity of labor supply  

rho = 0.075;   % time discount rate

% Production Function: F(k,n) = k^theta n^(1-theta)
theta = 0.3; % capital share

delta = 0.075; % depreciation rate

if delta == 0 
    disp(' ');
    disp(' *************************************************')
    disp('this will cause steady state investment to be zero')
    disp('so some quantities will not be well defined');
    disp(' *************************************************')
    disp('  ');
end
    

% labor and capital shares
sl = 1-theta ;
sk = theta ;
% elasticity of Fk with respecte to  k, e = -(k/Fk)(dFk/dk)
e = 1-theta ;


% Steady States
%**************

% capital output ratio
kovery = sk/(rho+delta);
% investment output ratio
xovery = delta*sk/(delta+rho);


% Variable Labor Supply case
% ***************************

% Elasticities of n(k,n):   vl(n)/vc(c) = Fn(k/n,1)
etak = theta/(theta+phi);
etac = sigma/(theta+phi);

% psi = (k/c)(dc/dk) elasticity of optimal consumption c(k) 
% Coefficient for 0= psi^2 a + psi b + c 

a = sl*etac + (rho+delta*sl)/(delta+rho)  ;
b = - ( 1 -(1-sl)*delta/(rho+delta) - sl + sl*etak )  - (e/sigma)*(1-sl)*etac ; 
c = (e/sigma)*(1-sl)*(etak-1);

psi1 = ( -b + sqrt( b^2 - 4*a*c)) / (2*a);
psi2 = ( -b - sqrt( b^2 - 4*a*c)) / (2*a);

psi = psi1;

% nu = (k/n)(dn/dk)  elasticity of optimal labor supply n(k)
nu = etak - etac*psi;


% mu = elasticity of output y(k,n(k))
% mu = (k/y)(dy/dk) = kFk/y + nFn/y ( (k/n)(dn/dk) + (c/n)(dn/dc) (k/c)(dc/dk) ) = sk + sl*( etak - etac*psi) 
mu = sk + sl*nu;

% alpha = (x/k)(dx/dk) elasticity of optimal investment x(k)
% c + x = y
% c/y psi + x/y alpha = mu or alpha = (mu - c/y psi )/(x/y)
alpha = (mu - (1-xovery)*psi)/xovery;

% (dkdot/dt) = lambda = speed of adjustment
lambda = (rho+delta)*(1 + nu*sl/(1-sl) - ( 1/(1-sl) - delta/(delta+rho) )*psi ) -delta;

% half life tstar= -log2/lambda ;
tstar = - log(2)/lambda;

% check of psi
check = (e/sigma)*(nu-1);
check = check / ( 1 + nu*sl/(1-sl) -delta/(delta+rho) - ( 1/(1-sl) - delta/(delta+rho) )*psi ) ;
check = psi - check;
checkzero = check;



% Fixed labor supply case
%************************

% psif = (k/c)(dc/dk) elasticity of optimal consumption c(k) 
% Coefficient for 0= psi^2 a + psi b + c 

a = (rho+delta*sl)/(delta+rho)  ;
b = - ( 1 -(1-sl)*delta/(rho+delta) - sl  )  ; 
c = -(e/sigma)*(1-sl);

psi1 = ( -b + sqrt( b^2 - 4*a*c)) / (2*a);
psi2 = ( -b - sqrt( b^2 - 4*a*c)) / (2*a);

psif = psi1;

% check of psif
checkf = -(e/sigma);
checkf = checkf / ( 1 - delta/(delta+rho) - ( 1/(1-sl) - delta/(delta+rho) )*psif ) ;
checkf = psif - checkf;
checkzerof = checkf;

% mu = elasticity of output y(k) ;
% mu = (k/y)(dy/dk) = kFk/y 
muf = sk ;

% alpha = (x/k)(dx/dk) elasticity of optimal investment x(k)
% c + x = y
% c/y psi + x/y alpha = mu or alpha = (mu - c/y psi )/(x/y)
alphaf = (muf - (1-xovery)*psif)/xovery;

% (dkdot/dt) = lambda = speed of adjustment
lambdaf = (rho+delta)*(1 - ( 1/(1-sl) - delta/(delta+rho) )*psif ) - delta;

% half life tstar= -log2/lambda ;
tstarf = - log(2)/lambdaf;


disp('parameters')
disp('**********')
disp(' ');
disp('sigma = CRRA coefficient')
disp(sigma)
disp(' ')
disp('phi = elasticity of v(c,1-n) w.r.t labor supply n ')
disp(phi);
disp(' ');
disp('rho = discount rate');
disp(rho);
disp(' ');
disp('theta = capital share')
disp(theta);
disp(' ');
disp('delta = depreciation rate')
disp(delta);
disp(' ');

disp('steady state quantities')
disp('***********************')
disp(' ')
disp('xovery = investment/output')
disp(xovery)
disp(' ')
disp('kovery = capital/output')
disp(kovery)
disp(' ')

disp('variable labor supply case')
disp('**************************')
disp(' ')
disp('next variable should be zero');
disp(checkzero);
disp(' ');
disp('psi  = (c/k)(dc/dk)');
disp(psi);
disp(' ' );
disp('nu  = (n/k)(dn/dk)');
disp(nu);
disp(' ' );
disp('mu  = (y/k)(dy/dk)');
disp(mu);
disp(' ' );
disp('alpha = (x/k)(dx/dk)');
disp(alpha);
disp(' ');
disp('lambda = dkdot/dk' );
disp(lambda);
disp(' ');
disp('tstar = half life' );
disp(tstar);
disp(' ');

disp('fixed labor supply case')
disp('***********************')
disp(' ')
disp('next variable should be zero');
disp(checkzerof);
disp(' ');
disp('psif  = (c/k)(dc/dk)');
disp(psif);
disp(' ' );
disp('muf  = (y/k)(dy/dk)');
disp(muf);
disp(' ' );
disp('alphaf = (x/k)(dx/dk)');
disp(alphaf);
disp(' ');
disp('lambdaf = dkdot/dk' );
disp(lambdaf);
disp(' ');
disp('tstarf = half life' );
disp(tstarf);
disp(' ')

%% Simulation

k0overkstar= 0.99 ; % initial capital / steady state capital
khat0 = log(k0overkstar); % initial value, as percentage of new ss 
% step is the size of time period, step = 0.25 is quarters if 1 is a year
step = 0.25;
capt = floor(tstar*3/step) ; % simulate up to 3 times half life
capt = max(capt,4);
% simulation for the variable labor supply caese

for s=1:capt
    t=(s-1)*step;
    time(s) = t;
    khat(s) = exp(lambda*t) * khat0 ;
    chat(s)= psi * khat(s);
    yhat(s) = mu * khat(s);
    xhat(s) = alpha * khat(s);
    nhat(s) = nu * khat(s); 
    what(s) = theta*(khat(s)-nhat(s));
end

figure(1);
plot(time,100*khat,time,100*chat,time,100*yhat,time,100*xhat,time,100*what,time,100*nhat);
title('Neoclassical Growth Model: adjustment to steady state, variable n')
legend('khat','chat','yhat','xhat','what','nhat')
ylabel('% deviation from steady state');
sxlabel = ['time periods in ', num2str(step), ' of a year'];
xlabel(sxlabel);


% simulation for the fixed labor supply case

for s=1:capt
    t=(s-1)*step;
    time(s) = t;
    khatf(s) = exp(lambdaf*t) * khat0 ;
    chatf(s)= psif * khatf(s);
    yhatf(s) = muf * khatf(s);
    xhatf(s) = alphaf * khatf(s);
    whatf(s) = theta*khatf(s);
end

figure(2);
plot(time,100*khatf,time,100*chatf,time,100*yhatf,time,100*xhatf,time,100*whatf);
title('Neoclassical Growth Model: adjustment to steady state, fixed n')
legend('khat','chat','yhat','xhat','what')
ylabel('% deviation from steady state');
sxlabel = ['time periods in ', num2str(step), ' of a year'];
xlabel(sxlabel);



% simulation for a permanent shock (!!!! ONLY VALID FOR SIGMA = 1 !!!)

if sigma ==1 ;
    
    inc = -log(k0overkstar);
    ktilde = khat + inc;
    ctilde = chat + inc;
    xtilde = xhat + inc;
    ytilde = yhat + inc;
    wtilde = what + inc;
    ntilde = nhat ;


    figure(3);
    plot(time,100*ktilde,time,100*ctilde,time,100*ytilde,time,100*xtilde,time,100*wtilde,time,100*ntilde);
    title('Neoclassical Growth Model: permanent shock to productivity, variable n')
    legend('ktilde','ctilde','ytilde','xtilde','wtilde','ntilde')
    ylabel('% deviation from initial values');
    sxlabel = ['time periods in ', num2str(step), ' of a year'];
    xlabel(sxlabel);
    
else
    disp('');
    disp('effect of permanent shock not computed because sigma not equal to 1');
    disp('*******************************************************************');
    disp(' ')
    
end
    

diary off;