#=
==========================================================================
% ACS_SS.M
% This program obtains the numerical values of the derivatives of the
% equilibrium conditions needed to compute the first order approximation
% for the New Keynesian model in ACS(2017).
%
% This function was created manually after calling num_eval. If the model
% changes this function has to be updated accordingly.
%
% Written by: Pablo Cuba-Borda, Federal Reserve Board
% Created: November 17, 2017. Washington D.C.
==========================================================================
=#

function get_order1_sgu(params_in)

phi    = params_in["phi"];
nu     = params_in["nu"];
tau    = params_in["tau"];
eta    = params_in["eta"];
psi1   = params_in["psi1"];
psi2   = params_in["psi2"];
beta   = params_in["beta"];
pi_ss  = params_in["pi_ss"];
gstar  = params_in["gstar"];
rho_z  = params_in["rho_z"];
rho_g  = params_in["rho_g"];
rho_r  = params_in["rho_r"];
rho_d  = params_in["rho_d"];

# Steady State for log-linearized variables
p1    = 0;
c1    = 0;
g1    = 0;
y1    = 0;
R1    = 0;
z1    = 0;
d1    = 0;
dy1   = 0;
m1    = 0;

# MAP STEADY STATE VALUES TO TIME t+1 VARIABLES
p2 = p1;
c2 = c1;
g2 = g1;
d2 = d1;
z2 = z1;
dy2 = dy1;
y2 = y1;
m2 = m1;

# MAP STEADY STATE VALUES TO TIME t-1 VARIABLES
y0 = y1;
R0 = R1;


nfx = [ 0      0         0               0                exp(R1 - d1 + d2 - p2 - z2 + c1*tau - c2*tau)  0;
 0      0         0                0  -beta*exp(d2 - d1 + dy2 + p2 + c1*tau - c2*tau)*(exp(p2) - 1)   0;
 0 -rho_r        0 psi2*(rho_r - 1)                                                            0 -1;
 0     0 exp(-g1)               0                                                            0  0;
 0     0        0               0                                                        -rho_d  0;
 0     0    -rho_g              0                                                            0  0;
 0     0        0           -rho_z                                                            0  0;
 1     0        0               0                                                            0  0;
 0     0        0               0                                                            0  0];


nfxp =[                                                0 -exp(R1 - d1 + d2 - p2 - z2 + c1*tau - c2*tau) 0 exp(R1 - d1 + d2 - p2 - z2 + c1*tau - c2*tau)              -exp(R1 - d1 + d2 - p2 - z2 + c1*tau - c2*tau) 0;
 -(exp(c1*tau + y1/eta)*(nu - 1))/(eta*nu*phi*pi_ss^2)                                              0 0                                             0 beta*exp(d2 - d1 + dy2 + p2 + c1*tau - c2*tau)*(exp(p2) - 1) 0;
                                                       0                                              1 0                                             0                                                           0 0;
                                           -exp(c1 - y1)                                              0 0                                             0                                                           0 0;
                                                       0                                              0 0                                             0                                                           1 0;
                                                       0                                              0 1                                             0                                                           0 0;
                                                       0                                              0 0                                             1                                                           0 0;
                                                      -1                                              0 0                                             0                                                           0 0;
                                                       0                                              0 0                                             0                                                           0 1];


nfy = [                                                               -tau*exp(R1 - d1 + d2 - p2 - z2 + c1*tau - c2*tau)               0                                                                                  0;
 beta*tau*exp(d2 - d1 + dy2 + p2 + c1*tau - c2*tau)*(exp(p2) - 1) - (tau*exp(c1*tau + y1/eta)*(nu - 1))/(nu*phi*pi_ss^2)               0 exp(p1)*(exp(p1)*(1/(2*nu) - 1) - 1/(2*nu)) + exp(p1)*(exp(p1) - 1)*(1/(2*nu) - 1);
                                                                                                                       0 psi2*(rho_r - 1)                                                                   psi1*(rho_r - 1);
                                                                                                            exp(c1 - y1)               0                                              gstar*phi*pi_ss^2*exp(p1)*(exp(p1) - 1);
                                                                                                                       0               0                                                                                  0;
                                                                                                                       0               0                                                                                  0;
                                                                                                                       0               0                                                                                  0;
                                                                                                                       0               1                                                                                  0;
                                                                                                                       0               0                                                                                  0];


nfyp =[         tau*exp(R1 - d1 + d2 - p2 - z2 + c1*tau - c2*tau)                                                            0                                                                       exp(R1 - d1 + d2 - p2 - z2 + c1*tau - c2*tau);
 -beta*tau*exp(d2 - d1 + dy2 + p2 + c1*tau - c2*tau)*(exp(p2) - 1)  beta*exp(d2 - d1 + dy2 + p2 + c1*tau - c2*tau)*(exp(p2) - 1) beta*exp(d2 - d1 + dy2 + p2 + c1*tau - c2*tau)*(exp(p2) - 1) + beta*exp(d2 - d1 + dy2 + p2 + c1*tau - c2*tau)*exp(p2);
                                                                0                                                            0                                                                                                                   0;
                                                                0                                                            0                                                                                                                   0;
                                                                0                                                            0                                                                                                                   0;
                                                                0                                                            0                                                                                                                   0;
                                                                0                                                            0                                                                                                                   0;
                                                                0                                                            0                                                                                                                   0;
                                                                0                                                            0                                                                                                                   0];

return nfx, nfxp, nfy, nfyp

end
