function F = ss_fun_GK(x)

%% Load model parameters
global M_
% NumberOfParameters = M_.param_nbr;                            % Number of deep parameters.
% for i = 1:NumberOfParameters                                  % Loop...
%   paramname = cell2mat(deblank(M_.param_names(i,:)));                   %    Get the name of parameter i. 
%   eval([ paramname ' = M_.params(' int2str(i) ');']);         %    Get the value of parameter i.
% end             



gy = M_.params(1);
alp = M_.params(2);
delta = M_.params(3);
zzetap = M_.params(4);
zzetaw = M_.params(5);
sigma_c = M_.params(6);
chi = M_.params(7);
gam = M_.params(8);
phiX = M_.params(9);
xip = M_.params(10);
xiw = M_.params(11);
gammap = M_.params(12);
gammaw = M_.params(13);
PIE_bar = M_.params(14);
g_bar = M_.params(15);
Rn_bar = M_.params(16);
rho_r = M_.params(17);
theta_pie = M_.params(18);
theta_y = M_.params(19);
theta_dy = M_.params(20);
rhoA = M_.params(21);
rhoG = M_.params(22);
rhoMCS = M_.params(23);
rhoMRSS = M_.params(24);
rhoMPS = M_.params(25);
rhoRPS = M_.params(26);
rhoIS = M_.params(27);
sigmaB = M_.params(28);
spreadtarg = M_.params(29);
levtarg = M_.params(30);
epsilon = M_.params(31);
xequitytarg = M_.params(32);
sigma_s = M_.params(33);
wr = M_.params(34);
apsi = M_.params(35);

H       = exp(x(1));
KYW     = exp(x(2));
thetaFF = exp(x(3));
xiB     = exp(x(4));

g       = g_bar/100;
betta   = (PIE_bar/100+1)/(Rn_bar/100+1)/(1+g)^(-sigma_c); % (273)

A       = 1;
MCS     = 1;
MRSS    = 1;
IS      = 1;
MPS     = 1;
RPS     = 1;

PIE     = PIE_bar/100+1;

S       = 0; % (368)
Sdash   = 0; % (369)
PIEw    = (1+g)*PIE; % (370)
Q       = 1; % (371)
LAMBDA  = betta*(1+g)^(-sigma_c); % (372)
PO      = ((1-xip*PIE^((gammap-1)*(1-zzetap)))/(1-xip))^(1/(1-zzetap)); % (373)
Deltap  = (1-xip)/(1-xip*PIE^(zzetap*(1-gammap)))*PO^(-zzetap); % (374)
MC      = (1-1/zzetap)*1/MCS*(1-xip*(1+g)*LAMBDA*PIE^(zzetap*(1-gammap)))/(1-xip*(1+g)*LAMBDA*PIE^((zzetap-1)*(1-gammap)))*PO; % (375)
PWP     = MC; % (376)
WO      = ((1-xiw*PIE^(gammaw*(1-zzetaw))*PIEw^(zzetaw-1))/(1-xiw))^(1/(1-zzetaw)); % (377)
Deltaw  = (1-xiw)/(1-xiw*PIEw^(zzetaw)/PIE^(gammaw*zzetaw))*WO^(-zzetaw); % (378)
Hd      = H/Deltaw; % (379)
rK      = (1-alp)*(1+g)/KYW*PWP; % (380)
YW      = A*Hd*(KYW/(1+g))^((1-alp)/alp); % (381)
K       = YW*KYW; % (382)
Y       = YW/Deltap; % (383)
I       = K/IS*(g+delta)/(1+g); % (384)
G       = gy*Y; % (385)
C       = Y-G-I; % (386)
JJw     = Hd/(1-xiw*LAMBDA*PIEw^zzetaw*PIE^(gammaw*(1-zzetaw)-1)); % (387)
W       = alp*PWP*YW/Hd; % (388)
Jw      = WO*W*JJw; % (389)
Wh      = (1-xiw*(1+g)*LAMBDA*PIEw^zzetaw/PIE^(gammaw*zzetaw))*(1-1/zzetaw)*Jw/(Hd*MRSS); % (390)
U       = (C*(1-chi/(1+g)))^(1-sigma_c)/(1-sigma_c)*exp((sigma_c-1)*H^(1+apsi)/(1+apsi)); % (391)
xequity = xequitytarg;
kappaFF = -epsilon/xequitytarg; % (410)
R       = 1/(LAMBDA); % (392)
RK      = rK+(1-delta); % (393)
spread  = RK-R; % (394)
Re      = 1/(LAMBDA); % (395)
q       = rK/(Re-1+delta); % (396)
Sasset  = K; % (397)
mue     = 0; % (398)
Nn      = xiB/(1+g)*(rK+(1-delta)*Q)*Sasset; % (399)
E       = xequity*Q*Sasset/q; % (400)
Thetax  = thetaFF*(1+epsilon*xequity+kappaFF*xequity^2/2); % (401)
No      = sigmaB*((rK+(1-delta-R)*Q)*Sasset-(rK+(1-delta-R)*q)*E+R*Nn)/(1+g-sigmaB*R); % (402)
N       = No+Nn; % (403)
phiB    = Q*Sasset/N; % (404)
omega   = (1-sigmaB)+sigmaB*Thetax*phiB; % (405)
nud     = LAMBDA*omega*R; % (406)
mus     = Thetax-nud/phiB-mue*xequity; % (407)

lambda  = (1-sigma_c)*U/(C-chi*C/(1+g))-betta*chi*(1-sigma_c)*U/(1+g)^(1-sigma_c)/((1+g)*C-chi*C); % (497)
UH      = -(1-sigma_c)*U*H^apsi;

F = [Wh+UH/lambda; mus-LAMBDA*omega*spread;  levtarg-Q*Sasset/(N+q*E); spreadtarg-spread];