% function [ys,check] = model_simple_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/>.

%UPD: previous to Dynare 4.6 M_ and options_ were called globally, but not anymore. Thus, we need to use these variables as inputs. 
function [ys,params, check] = simple_model_steadystate(ys,exo, M_, options_)

    % 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;
    
    options=optimset('display','off'); % set options for numerical solver

%% Enter model equations here
% endogenous variables (28 variables)
% some of the equationas are redundant in steady state, which one?? also,
    P = 1; % normalization??? THINK ABOUT IT
    Q = 1; %Equation 3): Price of capital
    pk = Q; % Equation 17: Definition of relevant capital price
    zeta= 0; % Equation 9): Nominal price setting FOC Non binding in steady state?
    P=1;
    lambda=1;
    eta=eta_bar;
    Y_0=0.5;
    z=0;z_f=0; v=0; v_f=0; 
    options=optimset('display','off'); 
    [Y,fval,exitflag]=fsolve(@(Y)find_y(Y, betta, tau, theta,eta_bar, delta, nu, alphaa, siggma, epsilon),Y_0,options);
    r=1/betta-1;
    r_taylor = r;
    R_gross=1+(1-tau)*r;
    k=alphaa/eta_bar*Y/((1-(1/R_gross-betta)*theta)/betta-1+delta);
    n=(Y/(k^alphaa))^(1/(1-nu));
    W=1/eta_bar*(1-nu)*Y/n;
    c=Y-k*delta;
    B=theta*k;
    d= Y-W*n -delta*k-B+B/R_gross;
    
    mu=1/R_gross-betta;

    V = d/(1-betta);  % Equation 10): Firm value
    G = G_bar;
    T = G + B*(1/R_gross - 1/(1+r)) ; % Equation 14): Government budget'
    y_f=Y;
    p_f = (Y/y_f)^((eta-1)/eta); % Equation 19): Demand equation for one firm

    Pi = Y^((eta-1)/eta)*y_f^(1/eta) - W/P*n; % Equation 21): Definition of earnings

    Qk = Q*k; % Equation  22.) Market value of capital
    invest=delta*k;

% -----------------------------------------------------------------------
% Variables for analysis
% -----------------------------------------------------------------------
delta_kn= 0;
delta_kn_current=0;
net_leverage= (B/P-Pi)/(Q*(1-delta)*k);
k_n=k/n;
k_n_current=k/n;
sales=P*y_f;
b_nominal        =B/(P*(1+r));


%% end own model equations
    for iter = 1:length(M_.params) %update parameters set in the file
        eval([ 'params(' num2str(iter) ',1) = ' 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(y, betta, tau, theta,eta_bar, delta, nu, alphaa, siggma, epsilon)
% compute y_f numerically

    r=1/betta-1;
    R_gross=1+(1-tau)*r;
 
    k = ((betta/eta_bar)*alphaa*y^((eta_bar-1)/eta_bar)*y^(1/eta_bar)) / (1-theta*(1/R_gross-betta)-betta*(1-delta));
    n=(y/(k^alphaa))^(1/(1-nu));
  
    w_d=1/eta_bar*(1-nu)*y/n;
    
    c=y-k*delta;
    
    w_s=(c)^(-siggma)/n^(1/epsilon);
    resid=w_s-w_d;

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

end
