%compute steady state of the model

%KR here M_ included in brackets
function [ys,params,check] = NAWMDynare_steadystate(ys,exo,M_,options_,print_dummy)

%KR commented:
%global M_

%num_para = M_.param_nbr; %size(M_.param_names,1);
%KR: here curly brackets changed
%for idxxx=1:1:num_para
%    paramname = M_.param_names{idxxx};
%    eval([ paramname ' = M_.params(' int2str(idxxx) ');']);
%end


for ii = 1:M_.param_nbr
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end

check = 0;

%steady state calcuations
Rbar              =nom_Rtarget;
Rgovbar           =GovRSpread/400+Rbar;
epsRPGbar         =Rgovbar/Rbar;
Rcorppbar         =CorpGovSpread/400+Rgovbar;
betta             =Pibar*gzbar^varsigma/Rgovbar;
epsRPKbar         =Rgovbar/Rcorppbar;
betastar          = betta;
sIM               = sIMC+sIMI;
sX                = sIM;
sC                = 1-sI-sG-sIG-sX+sIM;

%time to build - outlay share
b1=1-b0;


% computation of steady-state values for k, c, N and pI and calibration of psi, nuCbar, nuIbar and tauKbar
out_ss = NAWM_ss(betta, kappa, zetaa,...
    gzbar, alfa,  varphiHbar, varphiWbar,...
    muC, muI, sI, sG, sIMC, ...
    sIMI, sX, pIMbar, ...
    sIG,nuG,alphaG,...
    nuK,KG_K_target,alphaK,b0,KY_target,...
    dtx_ybar_target,sC,sTR,omega,...
    cJcISSratio,tauCrevYbar_target,tauWfrevYbar_target,...
    tauWhrevYbar_target,CorpGovSpread,Pibar,epsRPKbar,given_alphaG,given_alphaK,alphaG_target,alphaK_target,varsigma);


%solution
kbar              = out_ss(1);
cbar              = out_ss(2);
Nbar              = out_ss(3);
pIbar             = out_ss(4);
psi               = out_ss(5);
nuCbar            = out_ss(6);
nuIbar            = out_ss(7);
tauKbar           = out_ss(8);
kGbar             = out_ss(9);
deltaG            = out_ss(10);
delta             = out_ss(11);
tauNbar           = out_ss(12);
cJbar             = out_ss(13);
cIbar             = out_ss(14);
trJbar            = out_ss(15);
trIbar            = out_ss(16);
tauCbar           = out_ss(17);
tauWhbar          = out_ss(18);
tauWfbar          = out_ss(19);
alphaK            = out_ss(20);
alphaG            = out_ss(21);
CorpGovSpread     = out_ss(22);

%Interest rates etc
Rbar              =nom_Rtarget;
Rgovbar           =GovRSpread/400+Rbar;
epsRPGbar         =Rgovbar/Rbar;
Rcorppbar         =CorpGovSpread/400+Rgovbar;
betta             =Pibar*gzbar^varsigma/Rgovbar;
epsRPKbar         =Rgovbar/Rcorppbar;

% relative prices
pHbar             = ((1-(1-nuCbar)*pIMbar^(1-muC))/nuCbar)^(1/(1-muC));
pYbar             = pHbar;
pGbar             = pHbar;
pIGbar            = pHbar;
pXbar             = pYbar;

% real factor prices
rKbar             = 1/((1-tauKbar)*epsRPKbar)*(gzbar^varsigma/betta-(1-delta)*epsRPKbar-tauKbar*delta*epsRPKbar)*pIbar;
wbar              = (1-alfa)/alfa*kbar/Nbar/gzbar/(1+tauWfbar)*rKbar*(1+( (1-alphaK)/alphaK )^(1/nuK)*(kGbar/kbar)^( (nuK-1)/nuK ));

% final goods
qIbar             = (1-(1-delta)/gzbar)*kbar;
qCbar             = cbar;

% composite capital
if nuK>1-1e-4 && nuK<1+1e-4
    ktildebar     = alphaK^(-alphaK)*(1-alphaK)^(-(1-alphaK))*kbar^alphaK*kGbar^(1-alphaK); %cobb douglas
else
    ktildebar     =( alphaK^(1/nuK)*kbar^( (nuK-1)/nuK )+ (1-alphaK)^(1/nuK)*kGbar^( (nuK-1)/nuK )     )^( nuK/(nuK-1) );
end

%production
ybar              = gzbar^(-alfa)*ktildebar^alfa*Nbar^(1-alfa)-psi;

%intermediate goods
hCbar             = nuCbar*(pHbar/pCbar)^(-muC)*qCbar;
hIbar             = nuIbar*(pHbar/pIbar)^(-muI)*qIbar;
hIGbar            = sIG*(pYbar*ybar)/pIGbar;
hGbar             = sG*(pYbar*ybar)/pGbar;
hbar              = hCbar+hIbar+hGbar+hIGbar;
imCbar            = (1-nuCbar)*(pIMbar/pCbar)^(-muC)*qCbar;
imIbar            = (1-nuIbar)*(pIMbar/pIbar)^(-muI)*qIbar;


%expenditure components
cbar              = qCbar;
ibar              = qIbar;
gbar              = sG*(pYbar*ybar)/pGbar;
iGbar             = sIG*(pYbar*ybar)/pIGbar;
xbar              = ybar-hbar;
imbar             = imCbar+imIbar;

%fiscal variables
bbar              = sBbar*pYbar*ybar;
trbar             = sTR*pYbar*ybar;
dtxbar            = tauNbar*wbar*Nbar+tauKbar*kbar/gzbar*(rKbar-delta*pIbar);
tauKRevybar       = tauKbar*kbar/gzbar*(rKbar-delta*pIbar)/(pYbar*ybar);
tauNRevybar       = tauNbar*wbar*Nbar/(pYbar*ybar);
dtxbarwin         = dtxbar/(wbar*Nbar);
lumpbar           = (pGbar*gbar+pIGbar*iGbar+bbar/Pibar/gzbar+trbar-tauCbar*cbar-dtxbar-(tauWfbar+tauWhbar)*wbar*Nbar-bbar/Rgovbar);


%some useful checks for GDP ratios
lumpYbar=lumpbar/(pYbar*ybar);
trYbar=trbar/(pYbar*ybar);
gYbar=pGbar*gbar/(pYbar*ybar);
iGYbar=pIGbar*iGbar/(pYbar*ybar);
tauCrevYbar=tauCbar*cbar/(pYbar*ybar);
dtxrevYbar=dtxbar/(pYbar*ybar);
tauWfrevYbar=tauWfbar*wbar*Nbar/(pYbar*ybar);
tauWhrevYbar=tauWhbar*wbar*Nbar/(pYbar*ybar);
intpayYbar=(Rgovbar-1)/(Rgovbar*gzbar*Pibar)*sBbar;
intpayYbarAnnual=(Rgovbar^4-1)/(Rgovbar^4*gzbar^4*Pibar^4)*sBbar/4;
govexpYbar=trYbar+gYbar+iGYbar+intpayYbar;
govrevYbar=tauCrevYbar+dtxrevYbar+tauWfrevYbar+tauWhrevYbar;
govDefYbar=sBbar*((gzbar*Pibar-1)/(Pibar*gzbar*Rgovbar));
govDefYbarAnnual=sBbar*((gzbar^4*Pibar^4-1)/(Pibar^4*gzbar^4*Rgovbar^4))/4;

%total utility valued consumption
if nuG>1-1e-4 && nuG<1+1e-4
    cItildebar    = alphaG^(-alphaG)*(1-alphaG)^(-(1-alphaG))*cIbar^alphaG*gbar^(1-alphaG); %cobb douglas
else
    cItildebar   =( alphaG^(1/nuG)*cIbar^( (nuG-1)/nuG )+ (1-alphaG)^(1/nuG)*gbar^(  (nuG-1)/nuG))^(nuG/(nuG-1));
end

%abbreviations for marginal cost equation
nomin=(1-alphaK)^(1/nuK)*(kGbar/(pHbar/varphiHbar*gzbar*(ybar+psi)))^((nuK-1)/nuK);
denomin= ( ( (1-alfa)/( (1+tauWfbar)*wbar ))^( (alfa-1)*(nuK-1)/(alfa*nuK)) ) * ( (pHbar/varphiHbar)^(-(nuK-1)/(alfa*nuK)) );
xi1=nomin/denomin;
xi2=(nuK-1)/nuK*( ((1-alphaK)/alphaK )^(1/nuK)*(kGbar/kbar)^( (nuK-1)/nuK ))/(1+( (1-alphaK)/alphaK )^(1/nuK)*(kGbar/kbar)^( (nuK-1)/nuK));

%DTXWIN measurement equation
aabar=1+( (1-alphaK)/alphaK )^(1/nuK)*(kGbar/kbar)^( (nuK-1)/nuK);

%Capital to output ratios
kbYbar=kbar/(pYbar*ybar*4);
kGYbar=kGbar/(pYbar*ybar*4);

%non-Ricardian/Ricardian Consumption and Transfers
cJcIbar  = cJbar/cIbar;
trJtrIbar= trJbar/trIbar;

%interest rates
RNominalAnnual=400*(Rbar-1);
RgovNominalAnnual=400*(Rgovbar-1);
RKNominalPriorAgencyCostsAnnual=400*(Rcorppbar-1);
RKRealPriorAgencyCostsAnnual=400*(Rcorppbar/Pibar-1);
RgovRealAnnual=400*(Rgovbar/Pibar-1);
RKRealPostAgencyCostsAnnual=400*(Rcorppbar*epsRPKbar/Pibar-1);

%Foreign inflation and interest rate
PiYstarbar        = Pibar;
Rstarbar          = Rbar;
epsRPstarbar      =((1-epsRPGbar+gammaBstar)/gammaBstar)^gammaBstar;
RstarNominalAnnual=400*(Rstarbar-1);
RstarRealAnnual   =400*(((1-gammaBstar*epsRPstarbar^(1/gammaBstar)+gammaBstar))*(Rstarbar/Pibar))-400;

%coef. UIP
lamRPstar=(epsRPstarbar^(1/gammaBstar))/(1+gammaBstar-gammaBstar*epsRPstarbar^(1/gammaBstar));

%wage income as share of GDP
sWN=(wbar*Nbar)/(ybar*pYbar);

%set values since delta and deltaG are set
KG_K_target=kGbar/kbar;
KY_target=kbar/(pYbar*ybar)/4;

%convert AR(2) roots into coefs.
B1ystarystar=root1ystar+root2ystar;
B2ystarystar=-root1ystar*root2ystar;

check=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%constructing the vector ys for Dynare%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if check==0,%Dynare Switch indicating good steady state
    
    lambda=0;  Q=0;  rK=0; piC=0; s=0;   i=0;    w=0; N=0;    piH=0;     piX=0;    pY=0;   pIM=0; imI=0;     pI=0;  im=0;
    x=0; rgov=0;   GovIntPay=0;    r=0;  k=0;    y=0;  imC=0;   c=0;  cI=0;  cJ=0; trI=0;  trJ=0;
    bstar=0; pH=0; pX=0; piY=0;  piI=0; piIM=0;         epsilonRPG=0; epsilonRPK=0;
    epsilonRPstar=0; gz=0; epsilon=0; epsilonI=0; varphiW=0; varphiH=0;  varphiX=0; varphistar=0; epsilonIM=0;
    nustar=0; g=0; tr=0; tauC=0; tauWf=0; tauWh=0; iG=0; b=0;   tauN=0; lump=0; aiG=0;   tauK=0; piYstar=0; ystar=0; rstar=0; pO=0; pCX=0; dtxwin=0;
    cItilde=0; ktilde=0; kG=0; epsilonC=0; epsilonINP=0;
    
    Dy_obs=(gEbar-1+gzbar-1)*100; Dc_obs=(gEbar-1+gzbar-1)*100;
    Di_obs=(gEbar-1+gzbar-1)*100; Dx_obs=(gEbar-1+gzbar-1)*100; Dim_obs=(gEbar-1+gzbar-1)*100;
    E_obs=0; Dwn_obs=((gzbar-1)+(Pibar-1))*100; Dystar_obs=(gEbar-1+gzbar-1)*100;
    dtxwin_obs=0; E=0;
    piIM_obs=piIM+(Pibar-1)*100;
    piY_obs=piY+(Pibar-1)*100;
    piC_obs=piC+(Pibar-1)*100;
    r_obs=4*r+(Rbar-1)*400;
    piYstar_obs=piYstar+(PiYstarbar-1)*100;
    rstar_obs=4*rstar+(Rstarbar-1)*400;
    b_obs=0;
    
    %put steady states into ys vector

    %comment KR
    ys=NaN*ones(M_.orig_endo_nbr,1);
    
    
    %KR replace by curly brackets
    NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
        for ii = 1:NumberOfEndogenousVariables
             varname = M_.endo_names{ii};
            eval(['ys(' int2str(ii) ') = ' varname ';']);
        end


    
    
end;

% saving s.s. values that are used as params in static file
%KR replace by curly brackets

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



%print steady state if desired
if nargin==5,
    num_para=size(M_.params,1);
    disp(' ');
    disp('Parameters and Steady States:');
    disp('-----------------------------');
    disp(' ');
    for pppp=1:1:num_para
        disp([M_.param_names{pppp},'=',num2str(M_.params(pppp))]);
    end
end






