% 中等规模的DSGE @许文立，安徽大学经济学院/CIMERS
% 2019-2-24

function [ys,check]=dsge_steadystate(ys,exo)
% function [ys,check] = msacale_dsge_baseline_steadystate(ys,exo)
% 用数值算法来为msacale_dsge_baseline.mod计算稳态值
% 输入: 
%   - ys        [vector] 内生变量稳态值的初值向量
%   - exo       [vector] 外生变量值的向量
%
% 输出: 
%   - ys        [vector] 内生变量的稳态值向量
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impos restriction on parameters)
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;


% cdelta_1=1/cbeta-1+cxi_s*(1+cpi_target);

% 稳态表达式
% /*var mu    
% cgamma = 0.2828;
A=1;
Z=1;
v=1;
u=1;
pssi=cpsi_s;
tauc=ctau_c_s;
taun=ctau_n_s;
tauk=ctau_k_s;
pi=cpi_target;
varnothing = cvarnothing_s;
xi = cxi_s;
i=1/cbeta*(1+pi)-1;
% r_d = 0.0251256;
r_d= (1+pi)/cbeta  -1;
%negative
r_s = 1/(cbeta*(1+pi))-1;
mu = 1/xi*(cdelta_0+ cdelta_1 + 1/cbeta - 1);
q = mu + cbeta*(1+r_s);
% q = cdelta_0/xi;
% mu= q- cbeta*(1+r_s);
% R=(1+tauc)/(1-tauk)*(1/cbeta-(1-cdelta_0));
R= (q/(1-tauk)) * (1-cdelta_0+ mu*xi - 1/cbeta);

pi_star=(((1+pi)^(1-cepsilon_p)-cfi_p*(1+pi)^(czeta_p*(1-cepsilon_p)))/(1-cfi_p))^(1/(1-cepsilon_p))-1;
v_p=(1-cfi_p)*((1+pi_star)/(1+pi))^(-cepsilon_p)/(1-cfi_p*(1+pi)^(cepsilon_p*(1-czeta_p)));
mc=(cepsilon_p-1)/cepsilon_p*(1+pi_star)/(1+pi)*(1-cfi_p*cbeta*(1+pi)^(cepsilon_p*(1-czeta_p)))/(1-cfi_p*cbeta*(1+pi)^((1-cepsilon_p)*(czeta_p-1)));
w = (1-calpha)/calpha * R * (calpha * mc/R)^(1/(1-calpha));
w_star=w*((1-cfi_w*(1+pi)^((cepsilon_w-1)*(1-czeta_w)))/(1-cfi_w))^(1/(1-cepsilon_w));
N_d=((cepsilon_w-1)/cepsilon_w/cpsi_s*(1-taun)/(1+tauc)/(1-cb)*(1-cbeta*cb)*((1-comega)*(calpha*mc/R)^(calpha/(1-calpha))/v_p-cdelta_0*(calpha*mc/R)^(1/(1-calpha)))^(-1)*(w/w_star)^(-cepsilon_w*cchi)*((1-cfi_w*cbeta*(1+pi)^(cepsilon_w*(1+cchi)*(1-czeta_w)))/(1-cfi_w*cbeta*(1+pi)^((cepsilon_w-1)*(1-czeta_w))))*w_star)^(1/(1+cchi));
K_hat=(calpha*mc/R)^(1/(1-calpha))*N_d;
K=K_hat;
I=cdelta_0*K;
Y=A*K_hat^calpha*N_d^(1-calpha)/v_p;
G=comega*Y;
C=Y-G-I;
lamb=1/C*(1-cbeta*cb)/(1-cb)/(1+tauc);
% lamb=(1-cbeta*cb)/(C*(1-cb)*(1+tauc));
% d = cgamma* (1+pi)/lamb*1/(1+pi-cbeta*(1+r_d));
% d = 1/(1+pi);
d = (1+pi)/(cbeta * (1+r_d));
% lamb = mu/q; 
x1= lamb *mc*Y/(1-cfi_p*cbeta*(1+pi)^((1-czeta_p)*cepsilon_p));
x2= lamb *Y/(1-cfi_p*cbeta*(1+pi)^((czeta_p-1)*(1-cepsilon_p)));
h1=cpsi_s*(w/w_star)^(cepsilon_w*(1+cchi))*N_d^(1+cchi)/(1-cfi_w*cbeta*(1+pi)^(cepsilon_w*(1+cchi)*(1-czeta_w)));
h2=(1-taun)*lamb*(w/w_star)^cepsilon_w*N_d/(1-cfi_w*cbeta*(1+pi)^((cepsilon_w-1)*(1-czeta_w)));
M=I;
L=varnothing * M;
ita = M/(xi * K);
% L = xi * K +d;
% % cgamma = d*lamb*(1-cbeta*(1+r_d)/(1+pi))/v;
% q = 1/cdelta_1*(1-tauk)*R; 
%mu= q- cbeta*(1+r_s);
%R= q * cdelta_1/(1-tauk);
% mu=(1+tauc)*lamb;
% w=(1-calpha)*mc*(calpha*mc/R)^(calpha/(1-calpha));
% r_s = 1/cbeta-1
% mu = ((1/cbeta)-cdelta_1-(1-cdelta_0))/xi;
% mu=(1+tauc)*lamb;

check=0;


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
