/*
 * example1_dynamic.c : Computes dynamic model for Dynare
 *
 * Warning : this file is generated automatically by Dynare
 *           from model file (.mod)
 */
#include <math.h>
#include <stdlib.h>
#define max(a, b) (((a) > (b)) ? (a) : (b))
#define min(a, b) (((a) > (b)) ? (b) : (a))
double getPowerDeriv(double, double, int);
void Dynamic(double *y, double *x, int nb_row_x, double *params, double *steady_state, int it_, double *residual, double *g1, double *v2, double *v3)
{
  double lhs, rhs;

  /* Residual equations */
double T14 = pow(y[7],1+params[5]);
double T24 = exp(y[8]);
double T27 = exp(y[11]);
double T37 = params[2]*T27*y[9]+y[5]*(1-params[3]);
double T42 = exp(y[6]);
double T45 = T42*pow(y[0],params[2]);
double T46 = pow(y[7],1-params[2]);
double T47 = T45*T46;
double T88 = T42*getPowerDeriv(y[0],params[2],1);
double T90 = (-(T46*T88));
double T98 = getPowerDeriv(y[7],1+params[5],1);
double T100 = getPowerDeriv(y[7],1-params[2],1);
double T107 = (-(y[4]*T24*T27*y[10]))/(T27*y[10]*T27*y[10]);
double T112 = (-(params[0]*(T37*T107+y[4]*T24/(T27*y[10])*params[2]*T27*y[9])));
lhs =y[4]*params[4]*T14;
rhs =(1-params[2])*y[3];
residual[0]= lhs-rhs;
lhs =y[5];
rhs =params[0]*y[4]*T24/(T27*y[10])*T37;
residual[1]= lhs-rhs;
lhs =y[3];
rhs =T47;
residual[2]= lhs-rhs;
lhs =y[5];
rhs =T24*(y[3]-y[4])+(1-params[3])*y[0];
residual[3]= lhs-rhs;
lhs =y[6];
rhs =params[1]*y[1]+params[6]*y[2]+x[it_+0*nb_row_x];
residual[4]= lhs-rhs;
lhs =y[8];
rhs =y[1]*params[6]+params[1]*y[2]+x[it_+1*nb_row_x];
residual[5]= lhs-rhs;
  /* Jacobian  */
  if (g1 == NULL)
    return;
  else
    {
  g1[18]=(-(1-params[2]));
  g1[24]=params[4]*T14;
  g1[42]=y[4]*params[4]*T98;
  g1[55]=(-(params[0]*y[4]*T24/(T27*y[10])*params[2]*T27));
  g1[25]=(-(params[0]*T37*T24/(T27*y[10])));
  g1[61]=(-(params[0]*T37*(-(y[4]*T24*T27))/(T27*y[10]*T27*y[10])));
  g1[31]=1-params[0]*y[4]*T24/(T27*y[10])*(1-params[3]);
  g1[49]=(-(params[0]*y[4]*T24/(T27*y[10])*T37));
  g1[67]=T112;
  g1[20]=1;
  g1[2]=T90;
  g1[38]=(-T47);
  g1[44]=(-(T45*T100));
  g1[21]=(-T24);
  g1[27]=T24;
  g1[3]=(-(1-params[3]));
  g1[33]=1;
  g1[51]=(-(T24*(y[3]-y[4])));
  g1[10]=(-params[1]);
  g1[40]=1;
  g1[16]=(-params[6]);
  g1[76]=(-1);
  g1[11]=(-params[6]);
  g1[17]=(-params[1]);
  g1[53]=1;
  g1[83]=(-1);
    }
  /* Hessian for endogenous and exogenous variables */
  if (v2 == NULL)
    return;
  else
    {
  v2[0]=1;
  v2[48]=103;
  v2[96]=params[4]*T98;
  v2[1]=1;
  v2[49]=64;
  v2[97]=  v2[96];
  v2[2]=1;
  v2[50]=106;
  v2[98]=y[4]*params[4]*getPowerDeriv(y[7],1+params[5],2);
  v2[3]=2;
  v2[51]=66;
  v2[99]=(-(params[0]*params[2]*T27*T24/(T27*y[10])));
  v2[4]=2;
  v2[52]=131;
  v2[100]=  v2[99];
  v2[5]=2;
  v2[53]=150;
  v2[101]=(-(params[0]*params[2]*T27*(-(y[4]*T24*T27))/(T27*y[10]*T27*y[10])));
  v2[6]=2;
  v2[54]=137;
  v2[102]=  v2[101];
  v2[7]=2;
  v2[55]=145;
  v2[103]=(-(params[0]*T37*(-(T24*T27))/(T27*y[10]*T27*y[10])));
  v2[8]=2;
  v2[56]=67;
  v2[104]=  v2[103];
  v2[9]=2;
  v2[57]=151;
  v2[105]=(-(params[0]*T37*(-((-(y[4]*T24*T27))*(T27*T27*y[10]+T27*T27*y[10])))/(T27*y[10]*T27*y[10]*T27*y[10]*T27*y[10])));
  v2[10]=2;
  v2[58]=75;
  v2[106]=(-(params[0]*(1-params[3])*T24/(T27*y[10])));
  v2[11]=2;
  v2[59]=62;
  v2[107]=  v2[106];
  v2[12]=2;
  v2[60]=81;
  v2[108]=(-(params[0]*(1-params[3])*(-(y[4]*T24*T27))/(T27*y[10]*T27*y[10])));
  v2[13]=2;
  v2[61]=146;
  v2[109]=  v2[108];
  v2[14]=2;
  v2[62]=122;
  v2[110]=(-(params[0]*y[4]*T24/(T27*y[10])*params[2]*T27));
  v2[15]=2;
  v2[63]=135;
  v2[111]=  v2[110];
  v2[16]=2;
  v2[64]=117;
  v2[112]=(-(params[0]*T37*T24/(T27*y[10])));
  v2[17]=2;
  v2[65]=65;
  v2[113]=  v2[112];
  v2[18]=2;
  v2[66]=123;
  v2[114]=(-(params[0]*T37*(-(y[4]*T24*T27))/(T27*y[10]*T27*y[10])));
  v2[19]=2;
  v2[67]=149;
  v2[115]=  v2[114];
  v2[20]=2;
  v2[68]=118;
  v2[116]=(-(params[0]*y[4]*T24/(T27*y[10])*(1-params[3])));
  v2[21]=2;
  v2[69]=79;
  v2[117]=  v2[116];
  v2[22]=2;
  v2[70]=121;
  v2[118]=(-(params[0]*y[4]*T24/(T27*y[10])*T37));
  v2[23]=2;
  v2[71]=164;
  v2[119]=(-(params[0]*(y[4]*T24/(T27*y[10])*params[2]*T27+params[2]*T27*T107)));
  v2[24]=2;
  v2[72]=138;
  v2[120]=  v2[119];
  v2[25]=2;
  v2[73]=159;
  v2[121]=(-(params[0]*(T37*(-(T24*T27*y[10]))/(T27*y[10]*T27*y[10])+params[2]*T27*y[9]*T24/(T27*y[10]))));
  v2[26]=2;
  v2[74]=68;
  v2[122]=  v2[121];
  v2[27]=2;
  v2[75]=165;
  v2[123]=(-(params[0]*(T37*((-(y[4]*T24*T27))*T27*y[10]*T27*y[10]-(-(y[4]*T24*T27*y[10]))*(T27*T27*y[10]+T27*T27*y[10]))/(T27*y[10]*T27*y[10]*T27*y[10]*T27*y[10])+params[2]*T27*y[9]*(-(y[4]*T24*T27))/(T27*y[10]*T27*y[10]))));
  v2[28]=2;
  v2[76]=152;
  v2[124]=  v2[123];
  v2[29]=2;
  v2[77]=160;
  v2[125]=(-(params[0]*(1-params[3])*T107));
  v2[30]=2;
  v2[78]=82;
  v2[126]=  v2[125];
  v2[31]=2;
  v2[79]=163;
  v2[127]=T112;
  v2[32]=2;
  v2[80]=124;
  v2[128]=  v2[127];
  v2[33]=2;
  v2[81]=166;
  v2[129]=(-(params[0]*(params[2]*T27*y[9]*T107+T37*(T27*y[10]*T27*y[10]*(-(y[4]*T24*T27*y[10]))-(-(y[4]*T24*T27*y[10]))*(T27*y[10]*T27*y[10]+T27*y[10]*T27*y[10]))/(T27*y[10]*T27*y[10]*T27*y[10]*T27*y[10])+y[4]*T24/(T27*y[10])*params[2]*T27*y[9]+params[2]*T27*y[9]*T107)));
  v2[34]=3;
  v2[82]=1;
  v2[130]=(-(T46*T42*getPowerDeriv(y[0],params[2],2)));
  v2[35]=3;
  v2[83]=85;
  v2[131]=T90;
  v2[36]=3;
  v2[84]=7;
  v2[132]=  v2[131];
  v2[37]=3;
  v2[85]=91;
  v2[133]=(-T47);
  v2[38]=3;
  v2[86]=99;
  v2[134]=(-(T88*T100));
  v2[39]=3;
  v2[87]=8;
  v2[135]=  v2[134];
  v2[40]=3;
  v2[88]=105;
  v2[136]=(-(T45*T100));
  v2[41]=3;
  v2[89]=92;
  v2[137]=  v2[136];
  v2[42]=3;
  v2[90]=106;
  v2[138]=(-(T45*getPowerDeriv(y[7],1-params[2],2)));
  v2[43]=4;
  v2[91]=116;
  v2[139]=(-T24);
  v2[44]=4;
  v2[92]=51;
  v2[140]=  v2[139];
  v2[45]=4;
  v2[93]=117;
  v2[141]=T24;
  v2[46]=4;
  v2[94]=65;
  v2[142]=  v2[141];
  v2[47]=4;
  v2[95]=121;
  v2[143]=(-(T24*(y[3]-y[4])));
    }
}

/*
 * The k-th derivative of x^p
 */
double getPowerDeriv(double x, double p, int k)
{
#ifdef _MSC_VER
# define nearbyint(x) (fabs((x)-floor(x)) < fabs((x)-ceil(x)) ? floor(x) : ceil(x))
#endif
  if ( fabs(x) < 1e-012 && p > 0 && k > p && fabs(p-nearbyint(p)) < 1e-012 )
    return 0.0;
  else
    {
      int i = 0;
      double dxp = pow(x, p-k);
      for (; i<k; i++)
        dxp *= p--;
      return dxp;
    }
}