% Computation of the standard neoclassical growth model with linearization
% methods
%
% These codes are available on the World Wide Web at http://www.econ.upenn.edu/~jesusfv
% and http://www.econ.upenn.edu/~aruoba and http://www.econ.umn.edu/~rubio.
% They are part of the package prepared for the
% paper "Comparing Solution Methods for Dynamic Equilibrium Economies". 
% They may be freely reproduced for educational and research purposes, so long as they 
% are not altered, this copyright notice is reproduced with them, and they are not sold for profit.
% Consent of the corresponding author must be obtained before using all or any part of these codes
% in a publication. 
%
% Copyright © 2003 S. Boragan Aruoba, Jesus Fernandez-Villaverde and Juan
% F. Rubio-Ramirez

%----------------------------------------------------------------
% 0. Housekeeping
%----------------------------------------------------------------

clear all
format long
tic

%----------------------------------------------------------------
% 1. Calibration
%----------------------------------------------------------------

alpha  = 0.4;
beta   = 0.9896; 
delta  = 0.0196; 
theta  = 0.357; 
tau    = 2; 
rho    = 0.95; 
sigma  = 0.007; 
z = 0;

%----------------------------------------------------------------
% 2. Computation of the Steady State
%----------------------------------------------------------------

phi     = (1/alpha*(1/beta-1+delta))^(1/(1-alpha));
omega   = (phi^(1-alpha)-delta);
big_phi = theta/(1-theta)*(1-alpha)*phi^(-alpha);
  
k = big_phi/(omega+phi*big_phi);
l = phi*k;
c = omega*k;
r = alpha*(k^(alpha-1))*(l^(1-alpha));
y = k^alpha*(l^(1-alpha));
  
x = [k,l];

disp(' ')
disp('Steady State Values:')
disp(' ')
disp(sprintf(' y = %2.4f	k = %2.4f	l = %2.4f	c = %2.4f	r = %2.4f', y,k,l,c,r))

values_tau(1) = 2;
values_tau(2) = 10;
values_tau(3) = 50;
disp(' ')

for  different_tau = 1:3
  tau = values_tau(different_tau);
  
  disp(sprintf('tau =  %2.4f',tau))
  
  %----------------------------------------------------------------
  % 3. Building Matrices, linearization
  %----------------------------------------------------------------

  alpha1 = (theta*(1-tau)-1)/c;
  alpha2 = -(1-theta)*(1-tau)/(1-l);
  alpha3 = beta*((alpha*(1-alpha))/l)*k^(alpha-1)*l^(1-alpha)-(((1-tau)*(1-theta))/(1-l));
  alpha4 = alpha*beta*k^(alpha-1)*l^(1-alpha);
  alpha5 = beta*((alpha*(alpha-1))/k)*k^(alpha-1)*l^(1-alpha);
  alpha6 = -((alpha/l)+(1/((1-l))))*c;
  alpha7 = y*((1-alpha)/l);

  A = 1;
  B = (alpha/k)*c-y*(alpha/k)-(1-delta);
  C = alpha6-alpha7;
  D = c-y;
  F = 0;
  G = alpha1*alpha*c/k+alpha5;
  H = -alpha1*alpha*c/k;
  J = alpha1*alpha6+alpha3;
  K = -(alpha1*alpha6+alpha2);
  L = alpha1*c+alpha4;
  M = -alpha1*c;
  N = rho;

  %----------------------------------------------------------------
  % 4. UC
  %----------------------------------------------------------------

  P = 0.5*(-(B/A+K/J-(G*C)/(J*A))-sqrt((B/A+K/J-(G*C)/(J*A))^2-4*(((K*B-H*C)/(J*A)))));
  R = -(1/C)*(A*P+B);
  Q = (-D*(J*N+K)+C*L*N+C*M)/(A*J*N+A*K-C*G-C*J*R);
  S = (-A*L*N-A*M+D*G+D*J*R)/(A*J*N+A*K-C*G-C*J*R);
  disp(' ')
  disp('UC Approximation Policy Function for kp,l (k, z), levels:')
  disp(' ')
  disp([P, Q; R, S])
  Plineal(different_tau) = P;
  Qlineal(different_tau) = Q;
  Rlineal(different_tau) = R;
  Slineal(different_tau) = S;

  %---------------------------------------------------------------- 
  % 5. Building Matrices, loglinearization
  %----------------------------------------------------------------

  alpha1 = theta*(1-tau)-1;
  alpha2 = (1-theta)*(1-tau)*l/(1-l);
  alpha3 = alpha*beta*y/k;
  alpha4 = alpha3*(1-alpha);
  alpha5 = l/(1-l)+alpha;

  A = k;
  B = alpha*(c-y)-(1-delta)*k;
  C = (alpha-1)*y-alpha5*c;
  D = c-y;
  F = 0;
  G = alpha1*alpha-alpha4;
  H = -alpha*alpha1;
  J = alpha4-alpha1*alpha5-alpha2;
  K = alpha1*alpha5+alpha2;
  L = alpha3+alpha1;
  M = -alpha1;
  N = rho;

  %----------------------------------------------------------------
  % 6. UC
  %----------------------------------------------------------------

  P = 0.5*(-(B/A+K/J-(G*C)/(J*A))-sqrt((B/A+K/J-(G*C)/(J*A))^2-4*(((K*B-H*C)/(J*A)))));
  R = -(1/C)*(A*P+B);
  Q = (-D*(J*N+K)+C*L*N+C*M)/(A*J*N+A*K-C*G-C*J*R);
  S = (-A*L*N-A*M+D*G+D*J*R)/(A*J*N+A*K-C*G-C*J*R);
  disp(' ')
  disp('UC Approximation Policy Function for kp,l (k, z), loglinearization::')
  disp(' ')
  disp([P, Q; R, S])
  Plog(different_tau) = P;
  Qlog(different_tau) = Q;
  Rlog(different_tau) = R;
  Slog(different_tau) = S;
  
end

save policy_functions_linearization Plineal Qlineal Rlineal Slineal Plog Qlog Rlog Slog

toc