function [ys,params,check] = tank_steadystate(ys,exo,M_,options_)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2021 Alessandro Sardone.
% Permission is granted to copy, distribute and/or modify
% this document under the terms of the GNU Free
% Documentation License, Version 1.3 or any later version
% published by the Free Software Foundation; with no
% Invariant Sections, no Front-Cover Texts, and no
% Back-Cover Texts. A copy of the license is included in
% the section entitled GNU Free Documentation License".
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% computes the steady state for tank.mod and uses a numerical
% solver to do so
% Inputs: 
%   - ys        [vector] vector of ial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%   - M_        [structure] dynare model structure
%   - options_  [structure] dynare model options structure
%
%
% Output: 
%   - ys        [vector] vector of steady state values for the endogenous variables
%   - params    [vector] vector of parameter values of the model
%   - exo       [vector] vector of exogenous variables of the model

    check = false;
    options_ = options_; %#ok

    % read out parameters to access them with their name
    NumberOfParameters = M_.param_nbr;
    for ii = 1:NumberOfParameters
      paramname = char(M_.param_names(ii,:));
      eval([ paramname ' = M_.params(' int2str(ii) ');']);
      strpar.(paramname) = M_.params(ii);
    end
    
    % read out endogenous variables to access them with their name
    NumberOfEndo = M_.endo_nbr;
    strys.Init = nan;
    for ii = 1:NumberOfEndo
      varname = char(M_.endo_names(ii,:));
      strys.(varname) = ys(ii);
    end
    
    % read out exogenous variables to access them with their name
    NumberOfExo = M_.exo_nbr;
    strexo.Init = nan;
    for ii = 1:NumberOfExo
      exoname = char(M_.exo_names(ii,:));
      strexo.(exoname) = exo(ii);
    end

        %% end own model equations
        findGuesstemp = @(x)findGuess(x, strys, strpar, strexo);
        xstart = [1; 0.3; 0.5; 0.5];

        xopt = fsolve (@(x) findGuesstemp(x), xstart); 
        [~, strys, strpar, strexo] = findGuess(xopt, strys, strpar, strexo);
        [~, strys, strpar, strexo] = findGuess(xstart, strys, strpar, strexo);

 
    for iter = 1:length(M_.params) %update parameters set in the file
      M_.params(iter) = strpar.(char(M_.param_names(iter,:)));
    end
    params = M_.params;
    
    NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
    for ii = 1:NumberOfEndogenousVariables
      varname = char(M_.endo_names(ii,:));
      ys(ii) = strys.(varname);
    end
    
    NumberOfExogenousVariables = M_.exo_nbr; %auxiliary variables are set automatically
    for ii = 1:NumberOfExogenousVariables
      varname = char(M_.exo_names(ii,:));
      exo(ii) = strexo.(varname);
    end
end

function [feval, strys, strpar, strexo] = findGuess(x, strys, strpar, strexo)
    strys.y= x(1);   %output
    strys.E=x(2);    %total energy
    strys.pE=x(3);
    strys.em=x(4);

    %exogenous
    strys.lHM=1/3;                     % total labor HM
    strys.lS=1/3;                      % total labor S
    %strexo.muE=0;                     % abatement effort
    %strexo.muY=0;                     % abatement effort

    %steady state equations
    strys.q = 1;
    strys.r = strys.q/strpar.betta - (1-strpar.dellta);                     % real interest rate
    strys.Ay=1.00;                     % Y TFP
    strys.AE=1.00;                     % E TFP
    strys.emE = strpar.Qe*strys.em;
    strys.emY = strys.em - strys.emE;
    strpar.gammaE = strys.emE/(strys.E*(1-strexo.muE));
    strpar.gammaY = strys.emY/(strys.y(1-strexo.muY));
    strys.tauY=strpar.theta1*strpar.theta2*(strexo.muY^(strpar.theta2-1))/strpar.gammaY;                      % carbon tax
    strys.tauE=strpar.theta1*strpar.theta2*(strexo.muE^(strpar.theta2-1))/strpar.gammaE;
    strys.z=strys.y*strpar.theta1*strexo.muY^strpar.theta2 + strys.E*strpar.theta1*strexo.muE^strpar.theta2;
    strys.g = strys.tauE*strys.emE+strys.tauY*strys.emY;
    strys.mcY = (strpar.xxi -1)/strpar.xxi - (strpar.theta1*strexo.muY^strpar.theta2 + strys.tauY*strpar.gammaY*(1-strexo.muY)); % marginal cost lagrangian multiplier
    strys.kY = (strys.y* strpar.allpha * strys.mcY)/strys.r;
    strys.kE = strpar.allphaE * strys.E/strys.r * (strys.pE-strys.tauE*strpar.gammaE*(1-strexo.muE)-strpar.theta1*strexo.muE^strpar.theta2);
    strys.K = strys.kY + strys.kE;
    strys.inv = strpar.dellta*strys.K;
    strys.lE = (strys.E/(strys.AE*strys.kE^strpar.allphaE))^(1/(1-strpar.allphaE));
    strys.lY = strys.lHM*strpar.DELTAHH+strys.lS*(1-strpar.DELTAHH)-strys.lE;
    strys.w = (strys.y*strpar.uppsilon*strys.mcY)/strys.lY;
    strys.Ey = strys.y*(1-strpar.allpha-strpar.uppsilon)*strys.mcY/strys.pE;
    strys.Ec = strys.E*strpar.Qen;
    strys.c=strys.y-strys.z-strys.g-strys.inv;
    strys.pS = (1+strpar.omegaS*strys.pE^(strpar.eppsilon-1)-strpar.omegaS)^(1/(strpar.eppsilon-1));
    strys.pHM = (1+strpar.omegaHM*strys.pE^(strpar.eppsilon-1)-strpar.omegaHM)^(1/(strpar.eppsilon-1));
    strpar.gammaHM = strpar.omegaHM*(strys.pE/strys.pHM)^(strpar.eppsilon-1);
    strpar.gammaS = strpar.omegaS*(strys.pE/strys.pS)^(strpar.eppsilon-1);
    strys.xHM = strys.w*strys.lHM/strys.pHM;
    strpar.pssiHM =(strys.w*strys.xHM^(-strpar.phic))/(strys.lHM^strpar.phil*strys.pHM);
    strys.cHM = (1-strpar.gammaHM)*(1/strys.pHM)^(-strpar.eppsilon)*strys.xHM;
    strys.EcHM = strpar.gammaHM*(strys.pE/strys.pHM)^(-strpar.eppsilon)*strys.xHM;
    strys.EcS = (strys.Ec- strpar.DELTAHH * strys.EcHM)/(1-strpar.DELTAHH);
    strys.xS = strys.EcS/(strpar.gammaS) * (strys.pE/strys.pS)^strpar.eppsilon;
    strys.lambdaS = (strys.xS^(-strpar.phic))/strys.pS;
    strys.cS= (strys.c-strys.cHM*strpar.DELTAHH)/(1-strpar.DELTAHH);    
    strpar.pssiS = (strys.w*strys.lambdaS)/(strys.lS^strpar.phil);
    strys.x = strpar.DELTAHH*strys.xHM + (1-strpar.DELTAHH)*strys.xS;

   
    feval = nan(4,1);

    lhs = strys.y;
    rhs = strys.Ay*strys.Ey^(1-strpar.allpha-strpar.uppsilon)*strys.kY^strpar.allpha*strys.lY^strpar.uppsilon;
    feval(1) = 1 - lhs / rhs;
    lhs = strys.cS;
    rhs = (1-strpar.gammaS)*(1/strys.pS)^(-strpar.eppsilon)*strys.xS;
    feval(2) = 1 - lhs / rhs;
    lhs = strys.w;
    rhs = (strys.pE-strys.tauE*strpar.gammaE*(1-strexo.muE)-strpar.theta1*strexo.muE^strpar.theta2)*(1-strpar.allphaE)*strys.E/strys.lE;
    feval(3) = 1 - lhs / rhs;
    lhs = strys.Ec;
    rhs = strys.E - strys.Ey;
    feval(4) = 1 - lhs/rhs;
 

end