function [residual, g1, g2, g3] = fe_zlb_dynamic(y, x, params, steady_state, it_)
%
% Status : Computes dynamic model for Dynare
%
% Inputs :
%   y         [#dynamic variables by 1] double    vector of endogenous variables in the order stored
%                                                 in M_.lead_lag_incidence; see the Manual
%   x         [nperiods by M_.exo_nbr] double     matrix of exogenous variables (in declaration order)
%                                                 for all simulation periods
%   steady_state  [M_.endo_nbr by 1] double       vector of steady state values
%   params    [M_.param_nbr by 1] double          vector of parameter values in declaration order
%   it_       scalar double                       time period for exogenous variables for which to evaluate the model
%
% Outputs:
%   residual  [M_.endo_nbr by 1] double    vector of residuals of the dynamic model equations in order of 
%                                          declaration of the equations.
%                                          Dynare may prepend auxiliary equations, see M_.aux_vars
%   g1        [M_.endo_nbr by #dynamic variables] double    Jacobian matrix of the dynamic model equations;
%                                                           rows: equations in order of declaration
%                                                           columns: variables in order stored in M_.lead_lag_incidence followed by the ones in M_.exo_names
%   g2        [M_.endo_nbr by (#dynamic variables)^2] double   Hessian matrix of the dynamic model equations;
%                                                              rows: equations in order of declaration
%                                                              columns: variables in order stored in M_.lead_lag_incidence followed by the ones in M_.exo_names
%   g3        [M_.endo_nbr by (#dynamic variables)^3] double   Third order derivative matrix of the dynamic model equations;
%                                                              rows: equations in order of declaration
%                                                              columns: variables in order stored in M_.lead_lag_incidence followed by the ones in M_.exo_names
%
%
% Warning : this file is generated automatically by Dynare
%           from model file (.mod)

%
% Model equations
%

residual = zeros(15, 1);
lhs =y(5);
rhs =y(19)-params(2)*(y(13)-y(3)-y(7));
residual(1)= lhs-rhs;
lhs =y(3);
rhs =(1-params(16))*y(1)+params(16)*y(4);
residual(2)= lhs-rhs;
lhs =y(4);
rhs =y(3)*params(1)+y(5)*params(11)+y(8);
residual(3)= lhs-rhs;
residual(4) = y(13)+params(13)/400;
lhs =y(8);
rhs =x(it_, 1);
residual(5)= lhs-rhs;
lhs =y(7);
rhs =params(8)*y(2)+x(it_, 2);
residual(6)= lhs-rhs;
lhs =y(12);
rhs =y(21);
residual(7)= lhs-rhs;
lhs =y(6);
rhs =y(13)-y(3);
residual(8)= lhs-rhs;
lhs =y(11);
rhs =y(14)-y(21);
residual(9)= lhs-rhs;
lhs =y(9);
rhs =y(20)-params(2)*(y(14)-y(21)-y(7));
residual(10)= lhs-rhs;
lhs =y(10);
rhs =y(8)+params(1)*y(21)+params(11)*y(9);
residual(11)= lhs-rhs;
lhs =y(14);
rhs =(-params(13))/400;
residual(12)= lhs-rhs;
lhs =y(15);
rhs =y(18);
residual(13)= lhs-rhs;
lhs =y(16);
rhs =y(19);
residual(14)= lhs-rhs;
lhs =y(17);
rhs =y(4)-y(3);
residual(15)= lhs-rhs;
if nargout >= 2,
  g1 = zeros(15, 23);

  %
  % Jacobian matrix
  %

  g1(1,3)=(-params(2));
  g1(1,5)=1;
  g1(1,19)=(-1);
  g1(1,7)=(-params(2));
  g1(1,13)=params(2);
  g1(2,1)=(-(1-params(16)));
  g1(2,3)=1;
  g1(2,4)=(-params(16));
  g1(3,3)=(-params(1));
  g1(3,4)=1;
  g1(3,5)=(-params(11));
  g1(3,8)=(-1);
  g1(4,13)=1;
  g1(5,8)=1;
  g1(5,22)=(-1);
  g1(6,2)=(-params(8));
  g1(6,7)=1;
  g1(6,23)=(-1);
  g1(7,21)=(-1);
  g1(7,12)=1;
  g1(8,3)=1;
  g1(8,6)=1;
  g1(8,13)=(-1);
  g1(9,21)=1;
  g1(9,11)=1;
  g1(9,14)=(-1);
  g1(10,7)=(-params(2));
  g1(10,9)=1;
  g1(10,20)=(-1);
  g1(10,21)=(-params(2));
  g1(10,14)=params(2);
  g1(11,8)=(-1);
  g1(11,9)=(-params(11));
  g1(11,10)=1;
  g1(11,21)=(-params(1));
  g1(12,14)=1;
  g1(13,18)=(-1);
  g1(13,15)=1;
  g1(14,19)=(-1);
  g1(14,16)=1;
  g1(15,3)=1;
  g1(15,4)=(-1);
  g1(15,17)=1;

if nargout >= 3,
  %
  % Hessian matrix
  %

  g2 = sparse([],[],[],15,529);
if nargout >= 4,
  %
  % Third order derivatives
  %

  g3 = sparse([],[],[],15,12167);
end
end
end
end
