%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file includes the parameters and computes the steady state for the
% demand model.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% Housekeeping %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% specify parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

theta=0.8;
beta=0.96;
beta_hat=0.99;
delta=0; %0.05, 0.022
pi_ss=0;
xi_ssN=0.25;
eta=1;
epsilon=6;
lambda=0.5;
alpha=0.8;

phi_p=500; %500

gamma=0.9; % This is the inercia parameter for the collateral constraint

phi_h=0; % Adjustment costs of housing for the borrowers

rho_i=0.8; % Taylor rule parameter
phi_pi=1.5; % Taylor rule parameter (for now phi_x is zero, I won't even write it in the equations, this can be modified later but

rho_xi=0.9;

H_bar=2; %2

rho_xii=0.99; %0.99

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% Compute the Steady State %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% I call the parameters _ssN, the N to emphasize ('natural') that they are
% not in log, and we need to write them in logs.

n_ssN=(1/3);
n_hat_ssN=(1/3);
N_ssN=(1/3);

i_ssN=(1/beta_hat)-1;

Y_ssN=(N_ssN)^(alpha);

C_ssN=Y_ssN;

mult_ssN=1-beta*(1+i_ssN);

w_ssN=((epsilon-1)*(alpha/epsilon))/(Y_ssN^((1-alpha)/alpha));


one_pi_H_h_ssN=(w_ssN*n_ssN)/(( (1-mult_ssN*theta-beta*(1-delta))/(beta*xi_ssN) )+(1-theta)+(1+i_ssN)*theta-(1-delta));

c_ssN=w_ssN*n_ssN-(1+i_ssN)*theta*one_pi_H_h_ssN+(1-delta)*one_pi_H_h_ssN+theta*one_pi_H_h_ssN-one_pi_H_h_ssN;

c_hat_ssN=(C_ssN-lambda*c_ssN)/(1-lambda);

b_ssN=-theta*one_pi_H_h_ssN;

one_pi_H_h_hat_ssN=(beta_hat*xi_ssN*c_hat_ssN)/(1-beta_hat*(1-delta));


one_pi_H_Hbar=(1-lambda)*one_pi_H_h_hat_ssN+lambda*one_pi_H_h_ssN;

pi_H_ssN=(one_pi_H_Hbar/H_bar)-1;

h_ssN=one_pi_H_h_ssN/(1+pi_H_ssN);

h_hat_ssN=one_pi_H_h_hat_ssN/(1+pi_H_ssN);

% Back up psi and psi_hat
psi=(w_ssN/c_ssN)*(1/((n_ssN)^eta));
psi_hat=(w_ssN/c_hat_ssN)*(1/((n_hat_ssN)^eta));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% Report in logs the ones that are not rates %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n_ss=log(n_ssN);
n_hat_ss=log(n_hat_ssN);
N_ss=log(N_ssN);
i_ss=i_ssN;
Y_ss=log(Y_ssN);
w_ss=log(w_ssN);
C_ss=log(C_ssN);
h_ss=log(h_ssN);
c_ss=log(c_ssN);
c_hat_ss=log(c_hat_ssN);
h_hat_ss=log(h_hat_ssN);
mult_ss=log(mult_ssN);
pi_H_ss=pi_H_ssN;
%b_ssN is negative, so I'll compute it in levels (not in log)
b_ss=b_ssN;
xi_ss=xi_ssN;

denom_ss=(1/(1+100*rho_xii)); %This is just a rough value (obviously wrong), a rough guess.
% To be used just as a starting value for denom_ss, dynare should be able to find the right one in steady state
% (of course, note that it is invariant).

save param_nkpc theta beta beta_hat delta pi_ss pi_H_ss eta epsilon lambda alpha phi_p rho_i phi_pi ...
    n_ss n_hat_ss N_ss i_ss Y_ss w_ss C_ss h_ss c_ss c_hat_ss h_hat_ss psi psi_hat mult_ss phi_h ...
    xi_ss rho_xi H_bar b_ss rho_xii denom_ss gamma


dynare nkpc_v1

