function [ys,check] = NKfinal_one_stickiness_param_steadystate(ys,exo)
% function [ys,check] = Jermann_Quadrini_2012_NK_steadystate(ys,exo)
% computes the steady state for the Jermann_Quadrini_2012_NK.mod and uses a numerical
% solver to do so
% Inputs: 
%   - ys        [vector] vector of initial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%
% Output: 
%   - ys        [vector] vector of steady state values fpr the the endogenous variables
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impos restriction on parameters)
%
% Copyright (C) 2014-2016 Johannes Pfeifer
% 
%  This is free software: you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation, either version 3 of the License, or
%  (at your option) any later version.
% 
%  It is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
% 
%  For a copy of the GNU General Public License,
%  see <http://www.gnu.org/licenses/>.

global M_ 

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = deblank(M_.param_names(ii,:));
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end

% initialize indicator
check = 0;

%% Enter model equations here

% endogenous variables (41 variables)

% utilization rate
u_k  = 1;
u_pi = 1;

% two normalization is possible as Eq 6 and Eq 16 are redundant in SS.
P = 1; % normalization
Y = 1; % normalization 

ql = (betta*(cpn+gamma_lt)) / (1-betta*(1-gamma_lt)); % Equation 1 (HH bond: LT)
r  = 1/betta - 1; % Equation 1 (HH bond: ST)

Q_k = 1/(1-varrho*(1-iota)^2-2*varrho*(1-iota)*(1-betta)); 
Q_pi = 1/(1-varrho*(1-iota)^2-2*varrho*(1-iota)*(1-betta));% Equation 4 (Price of capital)
pk_k = Q_k; % Equation 18 (Price of capital in constraint)
pk_pi = Q_pi;
chat = 0; upsilonhat = 0; wopthat = 0; What = 0; % Equation 6 (SS deviation)

ql_k   = 1/(1+(1-tau_k)*(1/ql - 1)); ql_pi  = 1/(1+(1-tau_pi)*(1/ql - 1)); % Equation 19 (Gross price)

chi_k = 0; chi_pi = 0; % Equation 10 (Price FOC)

% Equation 9 (Bond demand)
mu_k = (ql_k - betta*((1-gamma_lt)*ql_k + (cpn + gamma_lt))) / ((1-netissue*betta*(1-gamma_lt))*ql);
mu_pi = (ql_pi - betta*((1-gamma_lt)*ql_pi + (cpn + gamma_lt))) / ((1-netissue*betta*(1-gamma_lt))*ql);

% shock variables (16 variables)

z = 0; z_k = 0; z_pi = 0;
zeta = 1; zeta_k = 1; zeta_pi = 1;
dd = 1; dd_k = 1; dd_pi = 1;
gamma = 1;
eta = eta_bar;
upsilon = upsilon_bar;
G = G_bar;
var_sigma = 0;
xi = xi_bar;
xi_k = xi_bar;
xi_pi = xi_bar;


% numerically find thek
tol = 1e-6;
diff = 1e6;
thek0 = 0.8;
iter = 0;
maxiter = 1000;

while diff > tol && iter < maxiter
    
    if ~isnan(thek) && iter == 0
        thek0 = thek;
    end

    % find y_k numerically

    if share_pi > 0 
        options=optimset('display','off'); % set options for numerical solver
        y_k_0 = 0.5;
        [y_k,fval,exitflag]=fsolve(@(y_k)find_y_k(y_k,betta,delta_k,delta_pi,eta,thek,thepi,theta_k,theta_pi,mu_k,mu_pi,xi,xi_k,xi_pi,share_pi,Y,Q_k,pk_k,Q_pi,ap1,a0,am1,am2,am3),y_k_0,options);

        k_k = ((betta/eta)*theta_k*Y^((eta-1)/eta)*y_k^(1/eta)) / (Q_k-(xi*xi_k)*mu_k*thek*(1-delta_k)*pk_k-Q_k*betta*(1-delta_k)); % Equation 3
        y_pi =  ((Y^(1/eta) - (1-share_pi)*y_k^(1/eta))/share_pi)^eta; % Equation 21
        k_pi = ((betta/eta)*(1+(xi*xi_pi)*mu_pi*thepi*((ap1/betta+a0+betta*am1+betta^2*am2+betta^3*am3)/(ap1+a0+am1+am2+am3))))*theta_pi*Y^((eta-1)/eta)*y_pi^(1/eta) / (Q_pi-Q_pi*betta*(1-delta_pi)); % Equation 3
        
        % Equation 17 (Output)
        n_k = (y_k / k_k^theta_k)^(1/(1-theta_k)); 
        n_pi = (y_pi / k_pi^theta_pi)^(1/(1-theta_pi)); 

        W = (1/eta) *Y^((eta-1)/eta)*(1-theta_k)*y_k^(1/eta) / n_k; % Equation 8 (Labor demand)    

    elseif share_pi == 0

        y_k  = Y;
        k_k = ((betta/eta)*theta_k*Y^((eta-1)/eta)*y_k^(1/eta)) / (Q_k-(xi*xi_k)*mu_k*thek*(1-delta_k)*pk_k-Q_k*betta*(1-delta_k));

        % Equation 17 (Output)
        n_k = (y_k / k_k^theta_k)^(1/(1-theta_k)); 

        W = (1/eta) *(1-theta_k)*Y^((eta-1)/eta)*y_k^(1/eta) / n_k; % Equation 8 (Labor demand)

        options=optimset('display','off'); % set options for numerical solver
        y_pi_0 = 0.5;
        [y_pi,fval,exitflag] = fsolve(@(y_pi)find_y_pi(y_pi,W,betta,delta_pi,eta,thepi,theta_pi,mu_pi,xi,xi_pi,Y,Q_pi,ap1,a0,am1,am2,am3),y_pi_0,options);

        k_pi = ((betta/eta)*(1+(xi*xi_pi)*mu_pi*thepi*((ap1/betta+a0+betta*am1+betta^2*am2+betta^3*am3)/(ap1+a0+am1+am2+am3))))*theta_pi*Y^((eta-1)/eta)*y_pi^(1/eta) / (Q_pi-Q_pi*betta*(1-delta_pi)); % Equation 3
        n_pi = (y_pi / k_pi^theta_pi)^(1/(1-theta_pi));    

    end

    invest_k = delta_k*k_k/(1-varrho*(1-iota)^2); invest_pi = delta_pi*k_pi/(1-varrho*(1-iota)^2); % Equation 5 (LOM capital)


    w_opt = W; % Equation 7 (Wage index)

    % Equation 12 (Enforcement)

    if netissue == 1
        b_pi = (((xi*xi_pi)*thepi * (Y^((eta-1)/eta)*y_pi^(1/eta) - W/P * n_pi)) / ql) / gamma_lt;
        if isnan(thek) || iter > 0
            thek = (((b_pi/y_pi)*y_k)/((xi*xi_k)*(1-delta_k)*k_k/ql)/gamma_lt);
        end
        b_k = (((xi*xi_k)*thek*(1-delta_k)*pk_k*k_k) / ql) / gamma_lt;    
    else
        b_pi = ((xi*xi_pi)*thepi * (Y^((eta-1)/eta)*y_pi^(1/eta) - W/P * n_pi)) / ql;
        if isnan(thek) || iter > 0
            thek = ((b_pi/y_pi)*y_k)/((xi*xi_k)*(1-delta_k)*k_k/ql);
        end
        b_k = ((xi*xi_k)*thek*(1-delta_k)*pk_k*k_k) / ql;
    end

    diff = abs((thek - thek0)/thek0);
    thek0 = thek;
    iter = iter + 1;

end

if iter == maxiter
   error('Cannot find thek numerically') 
end

% Equation 13 (Firm BC)
d_k = Y^((eta-1)/eta)*y_k^(1/eta) + ((gamma_lt * ql_k) - (cpn+gamma_lt))*b_k - W*n_k - invest_k;
d_pi = Y^((eta-1)/eta)*y_pi^(1/eta) + ((gamma_lt * ql_pi) - (cpn+gamma_lt))*b_pi - W*n_pi - invest_pi;

V_k = d_k/(1-betta); V_pi = d_pi/(1-betta); % Equation 11 (Firm value)

T = G + (1-share_pi)*gamma_lt*b_k*(ql_k - ql) + (share_pi)*gamma_lt*b_pi*(ql_pi - ql); % Equation 15 (Govt)

% Equation 21 (Agg.)
n = (1-share_pi)*n_k + share_pi*n_pi;
B = (1-share_pi)*b_k + share_pi*b_pi;
d = (1-share_pi)*d_k + share_pi*d_pi;
invest = (1-share_pi)*invest_k + share_pi*invest_pi;
B_s = 0;

c = W*n + ((cpn+gamma_lt) - gamma_lt*ql)*B + (1 - 1/(1+r))*B_s + d - T; % Equation 14 (HH BC)

P_k = (Y/y_k)^((eta-1)/eta); % Equation 20 (Demand equation)
if share_pi > 0
    P_pi = ((1-(1-share_pi)*P_k^(1/(1-eta)))/share_pi)^(1-eta);
elseif share_pi == 0
    P_pi = (Y/y_pi)^((eta-1)/eta);
end

vartheta_k = ((1/eta) * Y^((eta-1)/eta) * theta_k * y_k^(1/eta))/k_k;
vartheta_pi = (1+(xi*xi_pi)*mu_pi*thepi)*((1/eta) * Y^((eta-1)/eta) * theta_pi * y_pi^(1/eta))/k_pi;

% Variables for IRFs
c_pct           = log(c);

b_k_pct         = log((b_k - netissueirf*(1-gamma_lt)*b_k )/(P/ql) + 1e-6); % new debt issuance
b_pi_pct        = log((b_pi - netissueirf*(1-gamma_lt)*b_pi)/(P/ql) + 1e-6);
B_pct           = log((B - netissueirf*(1-gamma_lt)*B)/(P/ql) + 1e-6);

b_k_nom_pct     = log((b_k) + 1e-6);
b_pi_nom_pct    = log((b_pi) + 1e-6);
B_nom_pct       = log((B) + 1e-6);

invest_k_pct    = log(invest_k);
invest_pi_pct   = log(invest_pi);
invest_pct      = log((1-share_pi)*invest_k + share_pi*invest_pi);

y_k_pct         = log(y_k);
y_pi_pct        = log(y_pi);
Y_pct           = log(Y);

n_k_pct         = log(n_k);
n_pi_pct        = log(n_pi);
n_pct           = log(n);

pi_k_pct        = 0;
pi_pi_pct       = 0;
pi_pct           = 0;

z_k_pct = z_k;
z_pi_pct = z_pi;
z_pct = z;

zeta_k_pct = log(zeta_k);
zeta_pi_pct = log(zeta_pi);
zeta_pct = log(zeta);

dd_k_pct = log(dd_k);
dd_pi_pct = log(dd_pi);
dd_pct = log(dd);

eta_pct = log(eta);

mpn_k  = (1-theta_k)*y_k/n_k;
mpn_pi = (1-theta_pi)*y_pi/n_pi;

markup_k_pct    = log(P_k/(W/mpn_k));
markup_pi_pct   = log(P_pi/(W/mpn_pi));

Q_k_pct = log(Q_k);
Q_pi_pct = log(Q_pi);

Pi_k = Y^((eta-1)/eta)*y_k^(1/eta) - W/P*n_k;
Pi_pi = Y^((eta-1)/eta)*y_pi^(1/eta) - W/P*n_pi;
Pi = (1-share_pi)*Pi_k + share_pi*Pi_pi;

Pi_k_pct = log(Pi_k); Pi_pi_pct = log(Pi_pi); Pi_pct = log(Pi);

Qk_pi = Q_pi * k_pi;
V_Pi_pi = Pi_pi/(1-betta);
k_pi_rpl = k_pi;
SDF = betta;

pkk_k_pct = log(pk_k*k_k);
pkk_pi_pct = log(pk_pi*k_pi);
pkk_pct = log((1-share_pi)*pk_k*k_k + share_pi*pk_pi*k_pi);

W_pct = log(W);

eiss_k = -d_k; % equity issuance
eiss_pi = -d_pi;
eiss = -d;

% bissn_k = ((b_k - b_k)/P); % net debt issuance
% bissn_pi = ((b_pi - b_pi)/P);
% bissn = ((B - B)/P);

dissbyy_k = -d_k/(Y^((eta-1)/eta)*y_k^(1/eta));
dissbyy_pi = -d_pi/(Y^((eta-1)/eta)*y_pi^(1/eta));
dissbyy = -d/Y;

% bissbyy_k = ((b_k*ql - b_k)/P)/(Y^((eta-1)/eta)*y_k^(1/eta));
% bissbyy_pi = ((b_pi*ql - b_pi)/P)/(Y^((eta-1)/eta)*y_pi^(1/eta));
% bissbyy = ((B*ql - B)/P)/Y;

bissn_k = 0; % net debt issuance
bissn_pi = 0;
bissn = 0;

% observation equations
y_obs =0;
c_obs =0;
invest_obs =0;
pi_obs= 0;
r_obs= 0;
n_obs = 0;
W_obs= 0;
debt_obs = 0;
debt_k_obs = 0;
debt_pi_obs = 0;

%% end own model equations

for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'M_.params(' num2str(iter) ') = ' M_.param_names(iter,:) ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically

for ii = 1:NumberOfEndogenousVariables
  varname = deblank(M_.endo_names(ii,:));
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end

end

function resid = find_y_k(y_k,betta,delta_k,delta_pi,eta,thek,thepi,theta_k,theta_pi,mu_k,mu_pi,xi,xi_k,xi_pi,share_pi,Y,Q_k,pk_k,Q_pi,ap1,a0,am1,am2,am3)
% compute y_k numerically

k_k = ((betta/eta)*theta_k*Y^((eta-1)/eta)*y_k^(1/eta)) / (Q_k-(xi*xi_k)*mu_k*thek*(1-delta_k)*pk_k-Q_k*betta*(1-delta_k)); % Equation 3
y_pi =  ((Y^(1/eta) - (1-share_pi)*y_k^(1/eta))/share_pi)^eta; % Equation 21
k_pi = ((betta/eta)*(1+(xi*xi_pi)*mu_pi*thepi*((ap1/betta+a0+betta*am1+betta^2*am2+betta^3*am3)/(ap1+a0+am1+am2+am3))))*theta_pi*Y^((eta-1)/eta)*y_pi^(1/eta) / (Q_pi-Q_pi*betta*(1-delta_pi)); % Equation 3
n_k = (y_k / k_k^theta_k)^(1/(1-theta_k)); % Equation 17
n_pi = (y_pi / k_pi^theta_pi)^(1/(1-theta_pi)); % Equation 17

resid = ((1-theta_k) * y_k^(1/eta)/n_k - (1-theta_pi)*y_pi^(1/eta)/n_pi); % Equation 8

if imag(resid) ~= 0
    resid = 1e5;
end

end

function resid = find_y_pi(y_pi,W,betta,delta_pi,eta,thepi,theta_pi,mu_pi,xi,xi_pi,Y,Q_pi,ap1,a0,am1,am2,am3)
% compute y_pi numerically for share_pi = 0

k_pi = ((betta/eta)*(1+(xi*xi_pi)*mu_pi*thepi*((ap1/betta+a0+betta*am1+betta^2*am2+betta^3*am3)/(ap1+a0+am1+am2+am3))))*theta_pi*Y^((eta-1)/eta)*y_pi^(1/eta) / (Q_pi-Q_pi*betta*(1-delta_pi)); % Equation 3
n_pi = (y_pi / k_pi^theta_pi)^(1/(1-theta_pi));

resid = ((1/eta) *(1-theta_pi)*Y^((eta-1)/eta)*y_pi^(1/eta) / n_pi - W); % Equation 8 (labor demand)

if imag(resid) ~= 0
    resid = 1e5;
end

end



