function [ys,check] = rough_code_steadystate(ys,exo)

% Replicating CKM (steady-state file)
% 
% August 2018
% 
% This .m file is based on `NK_baseline_steadystate.m` from the dynare team...
%...github repository.
%
% This file computes the steady state of the model for CKM. It is necessary...
% ...because the steady state depends on the (endogenous) wedges, and therefore...
% ...must be computed numerically.

global M_ 

%Model-independent block
NumberOfParameters = M_.param_nbr; %Get number of params
for ii = 1:NumberOfParameters
  paramname = deblank(M_.param_names(ii,:)); %Read out parameter names
  eval([ paramname ' = M_.params(' int2str(ii) ');']); %Read out parameter values
end
% `check` is an indicator set to 0 if SS computation worked and 1 if not
check = 0;

%----------Model Equations----------
% Construct P0 matrix
P0=[P0_a; P0_n; P0_k; P0_g];

% Construct P matrix
P = [p_aa p_an p_ak p_ag
     p_na p_nn p_nk p_ng
     p_ka p_kn p_kk p_kg
     p_ga p_gn p_gk p_gg];
 
% Construct Sbar
Sbar=(eye(4,4)-P)\P0;

% Extract individual means
bar_tau_a   = Sbar(1);
bar_tau_n   = Sbar(2);
bar_tau_k   = Sbar(3);
bar_tau_g   = Sbar(4);

% 1 Mean of efficiency wedge
tau_a   = bar_tau_a;

% 2 Mean of labor wedge
tau_n   = bar_tau_n;

% 3 Mean of investment wedge
tau_k   = bar_tau_k;

% 4 Mean of government wedge
tau_g   = bar_tau_g;

% 5 Labor Share
n       = 1/3;                          
    
% 6 Interest rate definition
r       = (1 / bbeta) - (1 - ddelta);

% 7 Capital-output ratio (given wedges)
kkappa  = (alppha/r) * ((ggamma-1)/ggamma) * (1/(tau_a*tau_k));

% 8 Output (given K/Y and n)
y       = (tau_a^(ggamma/(ggamma-1))*(tau_k)^alppha*(tau_n)^(1-alppha)*kkappa^(alppha))^(1/(1-alppha)) * n;

% 9 Wage rate from firm optimization (given y, n, and wedges)
w       = (y/n)*(1-alppha)*((ggamma-1)/ggamma)*(1/(tau_a*tau_n));

% 10 k from firm optimization (given y, r, and wedges)
k       = y*(alppha/r)*((ggamma-1)/ggamma)*(1/(tau_a*tau_k));

% 11 Resource constraint defines c (given k and g wedge)
c       = y - (ddelta)*k - tau_g;

% 12 Investment (given k)
x       = k *ddelta;

%13 Define psi
ppsi    = w * (c^(-ssigma)) / (n^(1/pphi));

%k=log(k);
%w=log(w);
%c=log(c);
%y=log(y);
%x=log(x);
%l=log(l);

%----------End Model Equations----------

%Update parameters - based on dynare team code
for iter = 1:length(M_.params) %Update parameters
  eval([ 'M_.params(' num2str(iter) ') = ' M_.param_names(iter,:) ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %Set auxiliary vars automatically
for ii = 1:NumberOfEndogenousVariables
  varname = deblank(M_.endo_names(ii,:));
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end