%%%%%%%%%%%%% PINTUS-WEN 2011 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Replication file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% Housekeeping

clear all
clc


%% Calibration 

calibration = 'Calibration 2';

switch calibration
    case 'Calibration 1'
        
        L_fixed = 2;
        theta = 1;
        beta_tilde = 0.99;
        delta = 0.025;
        alpha = 0.35;
        gamma = 0.05;
        rho = 0.9;
        sigma_l = 1;
        sigma_w = 1;
        sigma_b = 4;
        beta    = 0.5;
        rho_a   = 0;
        
        %% Steady state

R_ss = 1/beta_tilde-1;
%%%%% compute numerically y_ss, k_ss and l_ss
param(1)=alpha;
param(2)=beta;
param(3)=delta;
param(4)=beta_tilde;
param(5)=sigma_w;
param(6)=sigma_l;
param(7)=gamma;
l_ss = 1; % I set land ratio equal to one; then I find a value for b consistent with this
param(8)=l_ss;
param(9)=L_fixed;

% pippo = [1 1 1];
% valori = steady_ykl(pippo,param);

% F(x)=0 where x = (y, k, l)
x0 = [1 1 1];
[x,val_eq] = fsolve(@(x) steady_ykl(x,param),x0);

% x0 = 0.5;
% [x,val_eq] = fsolve(@(x) steady_ykl(x,param),x0);

if abs(val_eq) > 0.0001
    error('steady state relations not satisfied!!!')
else
    disp('First check: STEADY STATE OK!!!!!')
end

b = x(1);
y_ss = (alpha*beta/(1-beta*(1-delta)))^(alpha/(1-alpha))*l_ss^(gamma/(1-alpha));
k_ss = (alpha*beta/(1-beta*(1-delta)))*y_ss;

clear x

r_ss = (1/beta_tilde)-1;
c_tilde_ss = beta*gamma*y_ss;
c_ss = y_ss-delta*k_ss-beta*gamma*y_ss;
q_ss = (1-beta_tilde)^(-1)*beta*gamma*(y_ss/l_ss); % ss price of land
b_ss = beta_tilde*q_ss*l_ss;
lambda_ss = (c_ss*(1-rho))^(-sigma_b);
lambda_tilde_ss = c_tilde_ss^(-sigma_l);
phi_ss = lambda_ss*(beta_tilde-beta);
l_tilde_ss = L_fixed-l_ss; 

disp('The steady state ratio of land allocated to lenders vs borrowers is:')
ratio_land = l_tilde_ss/l_ss
disp(' ')
disp('The steady state consumption level of lenders, relative to output, is:')
ratio_c_lender_GDP = c_tilde_ss/y_ss

% Check manually the steady state
eq16 = (c_ss-rho*c_ss)^(-sigma_b)-lambda_ss
eq17 = q_ss*lambda_ss-beta*q_ss*lambda_ss-beta*gamma*(y_ss/l_ss)*lambda_ss-theta*phi_ss*q_ss
eq18 = lambda_ss-(c_ss-rho*c_ss)^(-sigma_b)
eq19 = lambda_ss-beta*(1+r_ss)*lambda_ss-phi_ss*(1+r_ss)

eq13 = lambda_tilde_ss-c_tilde_ss^(-sigma_l)
eq14 = q_ss*lambda_tilde_ss-beta_tilde*q_ss*lambda_tilde_ss-beta_tilde*b*l_tilde_ss^(-sigma_w)
eq15 = lambda_tilde_ss-beta_tilde*(1+r_ss)*lambda_tilde_ss

somma = eq16+eq17+eq18+eq19+eq13+eq14+eq15;
if abs(somma)>0.0001
    error('steady state relations not satisfied!!!')
else
    disp('Second check: STEADY STATE OK!!!!!')
end
    

%% Call Dynare and solve the model


save param_pintus theta beta_tilde delta alpha gamma rho sigma_l sigma_w...
    sigma_b beta rho_a b y_ss k_ss l_ss c_tilde_ss c_ss q_ss b_ss lambda_ss...
    lambda_tilde_ss phi_ss l_tilde_ss L_fixed r_ss
        
dynare pintus_wen_dynarefile.mod noclearall  % running dynare

%%   Computing impulse response functions

% WITH DYNARE
% Impulse response to technology shock USING DYNARE
% Dynare will store IRFs for you under the following nomenclature: x_e denotes the impulse
% response of endogenous variable x to shock (varexo) e


% I adjust change_x = x_e/x_ss
% since I want percentage deviation,
% not absolute deviation from steady state
change_y        = (y_e_a)/y_ss;  % output
change_c_agg    = (c_agg_e_a)/c_ss; % aggregate consumption 
change_q        = (q_e_a)/q_ss;  % price of land
change_r        = (r_e_a);  % interest rate
change_l        = (l_e_a)/l_ss;  % borrower's land investment
change_k        = (k_e_a)/k_ss; % capital

T = length(change_y);
time = 1:1:T;
plot(time,change_y,'*-b',time,change_c_agg,'-k',time,change_q,'--r',time,change_r,':r',time,change_l,':k',time,change_k,'k')
title('Risk averse lender, \sigma =1 and i.i.d. TFP shock','fontsize',12,'fontweight','bold')
legend('Output','Agg. Cons.','Land Price','Interest rate','Land','Capital')

        
    case 'Calibration 2'
        
        L_fixed = 2;
        theta = 1;
        beta_tilde = 0.99;
        delta = 0.025;
        alpha = 0.35;
        gamma = 0.05;
        rho = 0.9;
        sigma_l = 1;
        sigma_w = 1;
        sigma_b = 2;
        rho_a   = 0.9;
        beta    = 0.8;
        
        %% Steady state

R_ss = 1/beta_tilde-1;
%%%%% compute numerically y_ss, k_ss and l_ss
param(1)=alpha;
param(2)=beta;
param(3)=delta;
param(4)=beta_tilde;
param(5)=sigma_w;
param(6)=sigma_l;
param(7)=gamma;
l_ss = 1; % I set land ratio equal to one; then I find a value for b consistent with this
param(8)=l_ss;
param(9)=L_fixed;

% pippo = [1 1 1];
% valori = steady_ykl(pippo,param);

% F(x)=0 where x = (y, k, l)
x0 = [1 1 1];
[x,val_eq] = fsolve(@(x) steady_ykl(x,param),x0);

if abs(val_eq) >0.0001
    error('steady state relations not satisfied!!!')
else
    disp('First check: STEADY STATE OK!!!!!')
end

b = x(1);
y_ss = (alpha*beta/(1-beta*(1-delta)))^(alpha/(1-alpha))*l_ss^(gamma/(1-alpha));
k_ss = (alpha*beta/(1-beta*(1-delta)))*y_ss;

clear x

r_ss = (1/beta_tilde)-1;
c_tilde_ss = beta*gamma*y_ss;
c_ss = y_ss-delta*k_ss-beta*gamma*y_ss;
q_ss = (1-beta_tilde)^(-1)*beta*gamma*(y_ss/l_ss); % ss price of land
b_ss = beta_tilde*q_ss*l_ss;
lambda_ss = (c_ss*(1-rho))^(-sigma_b);
lambda_tilde_ss = c_tilde_ss^(-sigma_l);
phi_ss = lambda_ss*(beta_tilde-beta);
l_tilde_ss = L_fixed-l_ss; 

disp('The steady state ratio of land allocated to lenders vs borrowers is:')
ratio_land = l_tilde_ss/l_ss
disp(' ')
disp('The steady state consumption level of lenders, relative to output, is:')
ratio_c_lender_GDP = c_tilde_ss/y_ss

% Check manually the steady state
eq16 = (c_ss-rho*c_ss)^(-sigma_b)-lambda_ss;
eq17 = q_ss*lambda_ss-beta*q_ss*lambda_ss-beta*gamma*(y_ss/l_ss)*lambda_ss-theta*phi_ss*q_ss;
eq18 = lambda_ss-(c_ss-rho*c_ss)^(-sigma_b);
eq19 = lambda_ss-beta*(1+r_ss)*lambda_ss-phi_ss*(1+r_ss);

eq13 = lambda_tilde_ss-c_tilde_ss^(-sigma_l);
eq14 = q_ss*lambda_tilde_ss-beta_tilde*q_ss*lambda_tilde_ss-beta_tilde*b*l_tilde_ss^(-sigma_w);
eq15 = lambda_tilde_ss-beta_tilde*(1+r_ss)*lambda_tilde_ss;

somma = eq16+eq17+eq18+eq19+eq13+eq14+eq15;
if abs(somma)>0.0001
    error('steady state relations not satisfied!!!')
else
    disp('Second check: STEADY STATE OK!!!!!')
end
    

%% Call Dynare and solve the model


save param_pintus theta beta_tilde delta alpha gamma rho sigma_l sigma_w...
    sigma_b beta rho_a b y_ss k_ss l_ss c_tilde_ss c_ss q_ss b_ss lambda_ss...
    lambda_tilde_ss phi_ss l_tilde_ss L_fixed r_ss
        
dynare pintus_wen_dynarefile.mod noclearall  % running dynare

%%   Computing impulse response functions

% WITH DYNARE
% Impulse response to technology shock USING DYNARE
% Dynare will store IRFs for you under the following nomenclature: x_e denotes the impulse
% response of endogenous variable x to shock (varexo) e


% I adjust change_x = x_e/x_ss
% since I want percentage deviation,
% not absolute deviation from steady state
change_y        = (y_e_a)/y_ss;  % output
change_c_agg    = (c_agg_e_a)/c_ss; % aggregate consumption 
change_q        = (q_e_a)/q_ss;  % price of land
change_r        = (r_e_a);  % interest rate
change_l        = (l_e_a)/l_ss;  % borrower's land investment
change_k        = (k_e_a)/k_ss; % capital


T = length(change_y);
time = 1:1:T;
plot(time,change_y,'*-b',time,change_c_agg,'-k',time,change_q,'--r',time,change_r,':r',time,change_l,':k',time,change_k,'k')
title('Risk averse lender, \sigma =1 and AR(1) TFP shock','fontsize',12,'fontweight','bold')
legend('Output','Agg. Cons.','Land Price','Interest rate','Land','Capital')

        
end