% 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_single_firm_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
%     betta = 0.99;
%     tau  = 0.1; 
%     alphaa = 0.25;
%     nu=0.75;
%     delta = 0.025;
%     theta=0.27;
%     siggma = 1.09;
%     epsilon = 0.44;
%     v=100;
%     phi = .25;
%     omega=0;
%     %omega = 0.158;
%     eta =1.11;
%     chi=12000;
%     rho_R = 0.785; 
%     nu_1 =2.5;
%     nu_2 = 0.125;
%     nu_3 = 0.0;
%     rho_xi = 0.503;
    pi=1;
    r_n=1/betta-1;
    R=1+(1-tau)*r_n;
    mu=1/R-betta;
    Q=1;
    M=eta;
    xi=0;
    z=1;

%guess of the x_c and x_u

%x_0=[3.85432, 1.17238];%for eta=1.3
x_0=[ 5.51152, 1.30631];%for eta=1.125
%options=optimset('display','off'); 
% options = optimset('TolFun',1e-30,'MaxFunEvals',1e5,'Maxiter',1e5,'Display','iter');
 options.TolX = 1e-45;
 options.MaxFunctionEvaluations=400;
[x,fval]=fsolve(@(x)find_phi_d_y_2(x, chi, eta, alphaa, betta, delta,nu, epsilon, theta, siggma, mu, Q, r_n, R,M,pi),x_0, options);
    k=x(1);
    y=x(2);
   
    psi_d=(betta*alphaa*y/M/k/(1-betta*(1-delta)-theta*mu)-1)/chi;
    psi_b=mu*(1+chi*psi_d)/chi;
    %k=betta*alphaa*y/M/(1-theta*chi*psi_b+chi*psi_d-betta*chi*psi_d+betta*delta*chi*psi_d);
 
    L= (y/(k^alphaa))^(1/nu);
    
    w= nu*y/M/L;
    
    c=(L^(1/epsilon)/w)^(-1/siggma);
    
    b=log(chi*psi_b)/chi+theta*k;
    
    T=b*(1/R- 1/(1+r_n));
    
    %d= y/M-w*L-delta*k-b+b/R-psi_b-psi_d;
    d= y/M-w*L-delta*k-b+b/R;
    invest=delta*k;
   
    p_s=betta/(1-betta)*d;

    r_r=r_n;
    dk=0;
   
%% 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_phi_d_y_2(x, chi, eta, alphaa, betta, delta,nu, epsilon, theta, siggma, mu, Q, r_n, R,M,pi)

    k=x(1);
    y=x(2);
    psi_d=(betta*alphaa*y/M/k/(1-betta*(1-delta)-theta*mu)-1)/chi;
    psi_b=mu*(1+chi*psi_d)/chi;
    
    L= (y/(k^alphaa))^(1/nu);
    
    w= nu*y/M/L;
    
    c_ls=(L^(1/epsilon)/w)^(-1/siggma);
    
    b=log(chi*psi_b)/chi+theta*k;
    
    T=b*(1/R-1/(1+r_n));
    
%     d_bc= y/M-w*L-delta*k-b+b/R;
%     c_bc= w*L +r_n/(1+r_n)*b+d_bc -psi_b-psi_d-T+(1-1/M)*y;
    d_bc= y/M-w*L-delta*k-b+b/R;

     c_bc= w*L +r_n/(1+r_n)*b+d_bc -psi_b-psi_d-T+(1-1/M)*y;

    d_const =-log(chi*psi_d)/chi;
    
    resid(1)=c_bc-c_ls;
    
        if imag(resid(1)) ~= 0
            resid(1) = 1e5;
        end
    
    resid(2)=d_bc-d_const;
        
        if imag(resid(2)) ~= 0
            resid(2) = 1e5;
        end


end


