function [residual, g1, g2, g3] = dynrbc_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(6, 1);
T37 = (-(exp(y(4))^(-params(4))));
T44 = exp(y(6))/exp(y(2));
T45 = T44-1;
T48 = 1+2*params(8)*T45/exp(y(2));
T53 = params(3)*exp(y(10))^(-params(4));
T61 = exp(y(6))^2;
T62 = (-exp(y(11)))/T61;
T69 = params(1)*exp(y(9))*exp(y(6))^(params(1)-1);
T70 = 1-params(2)-2*params(8)*(exp(y(11))/exp(y(6))-1)*T62+T69;
T91 = exp(y(3))*exp(y(2))^params(1);
lhs =T37*T48+T53*T70;
rhs =(-y(8))+params(3)*(1-params(2)+params(7))*y(12);
residual(1)= lhs-rhs;
lhs =exp(y(4))+exp(y(6))-exp(y(2))*(1-params(2))+params(8)*T45^2;
rhs =T91;
residual(2)= lhs-rhs;
lhs =exp(y(5));
rhs =exp(y(6))-exp(y(2))*(1-params(2));
residual(3)= lhs-rhs;
residual(4) = y(8);
lhs =y(3);
rhs =params(5)*y(1)+x(it_, 1);
residual(5)= lhs-rhs;
lhs =y(7);
rhs =y(2);
residual(6)= lhs-rhs;
if nargout >= 2,
  g1 = zeros(6, 13);

  %
  % Jacobian matrix
  %

T120 = (-(exp(y(6))*exp(y(2))))/(exp(y(2))*exp(y(2)));
  g1(1,9)=T53*T69;
  g1(1,4)=T48*(-(exp(y(4))*getPowerDeriv(exp(y(4)),(-params(4)),1)));
  g1(1,10)=T70*params(3)*exp(y(10))*getPowerDeriv(exp(y(10)),(-params(4)),1);
  g1(1,2)=T37*(exp(y(2))*2*params(8)*T120-exp(y(2))*2*params(8)*T45)/(exp(y(2))*exp(y(2)));
  g1(1,6)=T37*2*params(8)*T44/exp(y(2))+T53*((-(T62*2*params(8)*(-(exp(y(6))*exp(y(11))))/(exp(y(6))*exp(y(6)))+2*params(8)*(exp(y(11))/exp(y(6))-1)*(-((-exp(y(11)))*exp(y(6))*2*exp(y(6))))/(T61*T61)))+params(1)*exp(y(9))*exp(y(6))*getPowerDeriv(exp(y(6)),params(1)-1,1));
  g1(1,11)=T53*(-(2*params(8)*(exp(y(11))/exp(y(6))-1)*T62+T62*2*params(8)*exp(y(11))/exp(y(6))));
  g1(1,8)=1;
  g1(1,12)=(-(params(3)*(1-params(2)+params(7))));
  g1(2,3)=(-T91);
  g1(2,4)=exp(y(4));
  g1(2,2)=(-(exp(y(2))*(1-params(2))))+params(8)*T120*2*T45-exp(y(3))*exp(y(2))*getPowerDeriv(exp(y(2)),params(1),1);
  g1(2,6)=exp(y(6))+params(8)*T44*2*T45;
  g1(3,5)=exp(y(5));
  g1(3,2)=exp(y(2))*(1-params(2));
  g1(3,6)=(-exp(y(6)));
  g1(4,8)=1;
  g1(5,1)=(-params(5));
  g1(5,3)=1;
  g1(5,13)=(-1);
  g1(6,2)=(-1);
  g1(6,7)=1;

if nargout >= 3,
  %
  % Hessian matrix
  %

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

  g3 = sparse([],[],[],6,2197);
end
end
end
end
