Material for the Dynare Summer School 2018

deterministic-examples/rbc_basic.mod

+// Endogenous variables: consumption and capital
+var c k;
+
+// Exogenous variable: technology level
+varexo A;
+
+// Parameters declaration and calibration
+parameters alpha beta gamma delta;
+
+alpha=0.5;
+beta=0.95;
+gamma=0.5;
+delta=0.02;
+
+// Equilibrium conditions
+model;
+c + k = A*k(-1)^alpha + (1-delta)*k(-1); // Resource constraint
+c^(-gamma) = beta*c(+1)^(-gamma)*(alpha*A(+1)*k^(alpha-1) + 1 - delta); // Euler equation
+end;
+
+// Steady state (analytically solved)
+initval;
+A = 1;
+k = ((1-beta*(1-delta))/(beta*alpha*A))^(1/(alpha-1));
+c = A*k^alpha-delta*k;
+end;
+
+// Check that this is indeed the steady state
+steady;
+
+// Declare a positive technological shock in period 1
+shocks;
+var A;
+periods 1;
+values 1.2;
+end;
+
+// Prepare the deterministic simulation of the model over 100 periods
+perfect_foresight_setup(periods=100);
+
+// Perform the simulation
+perfect_foresight_solver;
+
+// Display the path of consumption
+rplot c;

deterministic-examples/rbc_det1.mod

+var k, y, L, c, A, a;
+
+varexo epsilon;
+
+parameters beta, theta, tau, alpha, psi, delta, rho, Astar;
+
+beta    =   0.9900;
+theta   =   0.3570;
+tau     =   2.0000;
+alpha   =   0.4500;
+psi     =  -0.1000;
+delta   =   0.0200;
+rho     =   0.8000;
+Astar   =   1.0000;
+
+model;
+ a = rho*a(-1) + epsilon;
+ A = Astar*exp(a);
+ y = A*(alpha*k(-1)^psi+(1-alpha)*L^psi)^(1/psi);
+ k = y-c + (1-delta)*k(-1);
+ (1-theta)/theta*c/(1-L) - (1-alpha)*(y/L)^(1-psi);
+ (c^theta*(1-L)^(1-theta))^(1-tau)/c = beta*(c(+1)^theta*(1-L(+1))^(1-theta))^(1-tau)/c(+1)*(alpha*(y(+1)/k)^(1-psi)+1-delta);
+end;
+
+steady_state_model;
+ a = epsilon/(1-rho);
+ A = Astar*exp(a);
+ Output_per_unit_of_Capital=((1/beta-1+delta)/alpha)^(1/(1-psi));
+ Consumption_per_unit_of_Capital=Output_per_unit_of_Capital-delta;
+ Labour_per_unit_of_Capital=(((Output_per_unit_of_Capital/A)^psi-alpha)/(1-alpha))^(1/psi);
+ Output_per_unit_of_Labour=Output_per_unit_of_Capital/Labour_per_unit_of_Capital;
+ Consumption_per_unit_of_Labour=Consumption_per_unit_of_Capital/Labour_per_unit_of_Capital;
+
+ % Compute steady state of the endogenous variables.
+ L=1/(1+Consumption_per_unit_of_Labour/((1-alpha)*theta/(1-theta)*Output_per_unit_of_Labour^(1-psi)));
+ c=Consumption_per_unit_of_Labour*L;
+ k=L/Labour_per_unit_of_Capital;
+ y=Output_per_unit_of_Capital*k;
+end;
+
+steady;
+
+ik = varlist_indices('k',M_.endo_names);
+kstar = oo_.steady_state(ik);
+
+histval;
+ k(0) = kstar/2;
+end;
+
+perfect_foresight_setup(periods=300);
+perfect_foresight_solver;
+
+rplot c;
+rplot k;

deterministic-examples/rbc_det2.mod

+var k, y, L, c, A, a;
+
+varexo epsilon;
+
+parameters beta, theta, tau, alpha, psi, delta, rho, Astar;
+
+beta    =   0.9900;
+theta   =   0.3570;
+tau     =   2.0000;
+alpha   =   0.4500;
+psi     =  -0.1000;
+delta   =   0.0200;
+rho     =   0.8000;
+Astar   =   1.0000;
+
+model;
+ a = rho*a(-1) + epsilon;
+ A = Astar*exp(a);
+ y = A*(alpha*k(-1)^psi+(1-alpha)*L^psi)^(1/psi);
+ k = y-c + (1-delta)*k(-1);
+ (1-theta)/theta*c/(1-L) - (1-alpha)*(y/L)^(1-psi);
+ (c^theta*(1-L)^(1-theta))^(1-tau)/c = beta*(c(+1)^theta*(1-L(+1))^(1-theta))^(1-tau)/c(+1)*(alpha*(y(+1)/k)^(1-psi)+1-delta);
+end;
+
+steady_state_model;
+ a = epsilon/(1-rho);
+ A = Astar*exp(a);
+ Output_per_unit_of_Capital=((1/beta-1+delta)/alpha)^(1/(1-psi));
+ Consumption_per_unit_of_Capital=Output_per_unit_of_Capital-delta;
+ Labour_per_unit_of_Capital=(((Output_per_unit_of_Capital/A)^psi-alpha)/(1-alpha))^(1/psi);
+ Output_per_unit_of_Labour=Output_per_unit_of_Capital/Labour_per_unit_of_Capital;
+ Consumption_per_unit_of_Labour=Consumption_per_unit_of_Capital/Labour_per_unit_of_Capital;
+
+ % Compute steady state of the endogenous variables.
+ L=1/(1+Consumption_per_unit_of_Labour/((1-alpha)*theta/(1-theta)*Output_per_unit_of_Labour^(1-psi)));
+ c=Consumption_per_unit_of_Labour*L;
+ k=L/Labour_per_unit_of_Capital;
+ y=Output_per_unit_of_Capital*k;
+end;
+
+steady;
+
+shocks;
+ var epsilon;
+ periods 1;
+ values -0.1;
+end;
+
+perfect_foresight_setup(periods=300);
+perfect_foresight_solver;
+
+rplot c;
+rplot k;

deterministic-examples/rbc_det3.mod

+var k, y, L, c, A, a;
+
+varexo epsilon;
+
+parameters beta, theta, tau, alpha, psi, delta, rho, Astar;
+
+beta    =   0.9900;
+theta   =   0.3570;
+tau     =   2.0000;
+alpha   =   0.4500;
+psi     =  -0.1000;
+delta   =   0.0200;
+rho     =   0.8000;
+Astar   =   1.0000;
+
+model;
+ a = rho*a(-1) + epsilon;
+ A = Astar*exp(a);
+ y = A*(alpha*k(-1)^psi+(1-alpha)*L^psi)^(1/psi);
+ k = y-c + (1-delta)*k(-1);
+ (1-theta)/theta*c/(1-L) - (1-alpha)*(y/L)^(1-psi);
+ (c^theta*(1-L)^(1-theta))^(1-tau)/c = beta*(c(+1)^theta*(1-L(+1))^(1-theta))^(1-tau)/c(+1)*(alpha*(y(+1)/k)^(1-psi)+1-delta);
+end;
+
+steady_state_model;
+ a = epsilon/(1-rho);
+ A = Astar*exp(a);
+ Output_per_unit_of_Capital=((1/beta-1+delta)/alpha)^(1/(1-psi));
+ Consumption_per_unit_of_Capital=Output_per_unit_of_Capital-delta;
+ Labour_per_unit_of_Capital=(((Output_per_unit_of_Capital/A)^psi-alpha)/(1-alpha))^(1/psi);
+ Output_per_unit_of_Labour=Output_per_unit_of_Capital/Labour_per_unit_of_Capital;
+ Consumption_per_unit_of_Labour=Consumption_per_unit_of_Capital/Labour_per_unit_of_Capital;
+
+ % Compute steady state of the endogenous variables.
+ L=1/(1+Consumption_per_unit_of_Labour/((1-alpha)*theta/(1-theta)*Output_per_unit_of_Labour^(1-psi)));
+ c=Consumption_per_unit_of_Labour*L;
+ k=L/Labour_per_unit_of_Capital;
+ y=Output_per_unit_of_Capital*k;
+end;
+
+steady;
+
+shocks;
+ var epsilon;
+ periods 4, 5:8;
+ values 0.04, 0.01;
+end;
+
+perfect_foresight_setup(periods=300);
+perfect_foresight_solver;
+
+rplot c;
+rplot k;

deterministic-examples/rbc_det4.mod

+var k, y, L, c, A, a;
+
+varexo epsilon;
+
+parameters beta, theta, tau, alpha, psi, delta, rho, Astar;
+
+beta    =   0.9900;
+theta   =   0.3570;
+tau     =   2.0000;
+alpha   =   0.4500;
+psi     =  -0.1000;
+delta   =   0.0200;
+rho     =   0.8000;
+Astar   =   1.0000;
+
+model;
+ a = rho*a(-1) + epsilon;
+ A = Astar*exp(a);
+ y = A*(alpha*k(-1)^psi+(1-alpha)*L^psi)^(1/psi);
+ k = y-c + (1-delta)*k(-1);
+ (1-theta)/theta*c/(1-L) - (1-alpha)*(y/L)^(1-psi);
+ (c^theta*(1-L)^(1-theta))^(1-tau)/c = beta*(c(+1)^theta*(1-L(+1))^(1-theta))^(1-tau)/c(+1)*(alpha*(y(+1)/k)^(1-psi)+1-delta);
+end;
+
+steady_state_model;
+ a = epsilon/(1-rho);
+ A = Astar*exp(a);
+ Output_per_unit_of_Capital=((1/beta-1+delta)/alpha)^(1/(1-psi));
+ Consumption_per_unit_of_Capital=Output_per_unit_of_Capital-delta;
+ Labour_per_unit_of_Capital=(((Output_per_unit_of_Capital/A)^psi-alpha)/(1-alpha))^(1/psi);
+ Output_per_unit_of_Labour=Output_per_unit_of_Capital/Labour_per_unit_of_Capital;
+ Consumption_per_unit_of_Labour=Consumption_per_unit_of_Capital/Labour_per_unit_of_Capital;
+
+ % Compute steady state of the endogenous variables.
+ L=1/(1+Consumption_per_unit_of_Labour/((1-alpha)*theta/(1-theta)*Output_per_unit_of_Labour^(1-psi)));
+ c=Consumption_per_unit_of_Labour*L;
+ k=L/Labour_per_unit_of_Capital;
+ y=Output_per_unit_of_Capital*k;
+end;
+
+initval;
+ epsilon = 0;
+end;
+
+steady;
+
+endval;
+ epsilon = (1-rho)*log(1.05);
+end;
+
+steady;
+
+perfect_foresight_setup(periods=300);
+perfect_foresight_solver;
+
+rplot c;
+rplot k;

deterministic-examples/rbc_det5.mod

+var k, y, L, c, A, a;
+
+varexo epsilon;
+
+parameters beta, theta, tau, alpha, psi, delta, rho, Astar;
+
+beta    =   0.9900;
+theta   =   0.3570;
+tau     =   2.0000;
+alpha   =   0.4500;
+psi     =  -0.1000;
+delta   =   0.0200;
+rho     =   0.8000;
+Astar   =   1.0000;
+
+model;
+ a = rho*a(-1) + epsilon;
+ A = Astar*exp(a);
+ y = A*(alpha*k(-1)^psi+(1-alpha)*L^psi)^(1/psi);
+ k = y-c + (1-delta)*k(-1);
+ (1-theta)/theta*c/(1-L) - (1-alpha)*(y/L)^(1-psi);
+ (c^theta*(1-L)^(1-theta))^(1-tau)/c = beta*(c(+1)^theta*(1-L(+1))^(1-theta))^(1-tau)/c(+1)*(alpha*(y(+1)/k)^(1-psi)+1-delta);
+end;
+
+steady_state_model;
+ a = epsilon/(1-rho);
+ A = Astar*exp(a);
+ Output_per_unit_of_Capital=((1/beta-1+delta)/alpha)^(1/(1-psi));
+ Consumption_per_unit_of_Capital=Output_per_unit_of_Capital-delta;
+ Labour_per_unit_of_Capital=(((Output_per_unit_of_Capital/A)^psi-alpha)/(1-alpha))^(1/psi);
+ Output_per_unit_of_Labour=Output_per_unit_of_Capital/Labour_per_unit_of_Capital;
+ Consumption_per_unit_of_Labour=Consumption_per_unit_of_Capital/Labour_per_unit_of_Capital;
+
+ % Compute steady state of the endogenous variables.
+ L=1/(1+Consumption_per_unit_of_Labour/((1-alpha)*theta/(1-theta)*Output_per_unit_of_Labour^(1-psi)));
+ c=Consumption_per_unit_of_Labour*L;
+ k=L/Labour_per_unit_of_Capital;
+ y=Output_per_unit_of_Capital*k;
+end;
+
+initval;
+ epsilon = 0;
+end;
+
+steady;
+
+endval;
+ epsilon = (1-rho)*log(1.05);
+end;
+
+steady;
+
+shocks;
+ var epsilon;
+ periods 1:5;
+ values 0;
+end;
+
+perfect_foresight_setup(periods=300);
+perfect_foresight_solver;
+
+rplot c;
+rplot k;

deterministic-examples/rbc_ep.mod

+var k, y, L, c, A, a;
+
+varexo epsilon;
+
+parameters beta, theta, tau, alpha, psi, delta, rho, Astar;
+
+beta    =   0.9900;
+theta   =   0.3570;
+tau     =   2.0000;
+alpha   =   0.4500;
+psi     =  -0.1000;
+delta   =   0.0200;
+rho     =   0.8000;
+Astar   =   1.0000;
+
+model;
+ a = rho*a(-1) + epsilon;
+ A = Astar*exp(a);
+ y = A*(alpha*k(-1)^psi+(1-alpha)*L^psi)^(1/psi);
+ k = y-c + (1-delta)*k(-1);
+ (1-theta)/theta*c/(1-L) - (1-alpha)*(y/L)^(1-psi);
+ (c^theta*(1-L)^(1-theta))^(1-tau)/c = beta*(c(+1)^theta*(1-L(+1))^(1-theta))^(1-tau)/c(+1)*(alpha*(y(+1)/k)^(1-psi)+1-delta);
+end;
+
+steady_state_model;
+ a = epsilon/(1-rho);
+ A = Astar*exp(a);
+ Output_per_unit_of_Capital=((1/beta-1+delta)/alpha)^(1/(1-psi));
+ Consumption_per_unit_of_Capital=Output_per_unit_of_Capital-delta;
+ Labour_per_unit_of_Capital=(((Output_per_unit_of_Capital/A)^psi-alpha)/(1-alpha))^(1/psi);
+ Output_per_unit_of_Labour=Output_per_unit_of_Capital/Labour_per_unit_of_Capital;
+ Consumption_per_unit_of_Labour=Consumption_per_unit_of_Capital/Labour_per_unit_of_Capital;
+
+ % Compute steady state of the endogenous variables.
+ L=1/(1+Consumption_per_unit_of_Labour/((1-alpha)*theta/(1-theta)*Output_per_unit_of_Labour^(1-psi)));
+ c=Consumption_per_unit_of_Labour*L;
+ k=L/Labour_per_unit_of_Capital;
+ y=Output_per_unit_of_Capital*k;
+end;
+
+steady;
+
+// Declare shocks as in a stochastic setup
+shocks;
+ var epsilon;
+ stderr 0.02;
+end;
+
+extended_path(periods=300);
+
+// Plot 20 first periods of consumption
+ic = varlist_indices('c',M_.endo_names);
+plot(oo_.endo_simul(ic, 1:21));

deterministic-examples/rbc_unexpected.mod

+var k, y, L, c, A, a;
+
+varexo epsilon;
+
+parameters beta, theta, tau, alpha, psi, delta, rho, Astar;
+
+beta    =   0.9900;
+theta   =   0.3570;
+tau     =   2.0000;
+alpha   =   0.4500;
+psi     =  -0.1000;
+delta   =   0.0200;
+rho     =   0.8000;
+Astar   =   1.0000;
+
+model;
+ a = rho*a(-1) + epsilon;
+ A = Astar*exp(a);
+ y = A*(alpha*k(-1)^psi+(1-alpha)*L^psi)^(1/psi);
+ k = y-c + (1-delta)*k(-1);
+ (1-theta)/theta*c/(1-L) - (1-alpha)*(y/L)^(1-psi);
+ (c^theta*(1-L)^(1-theta))^(1-tau)/c = beta*(c(+1)^theta*(1-L(+1))^(1-theta))^(1-tau)/c(+1)*(alpha*(y(+1)/k)^(1-psi)+1-delta);
+end;
+
+steady_state_model;
+ a = epsilon/(1-rho);
+ A = Astar*exp(a);
+ Output_per_unit_of_Capital=((1/beta-1+delta)/alpha)^(1/(1-psi));
+ Consumption_per_unit_of_Capital=Output_per_unit_of_Capital-delta;
+ Labour_per_unit_of_Capital=(((Output_per_unit_of_Capital/A)^psi-alpha)/(1-alpha))^(1/psi);
+ Output_per_unit_of_Labour=Output_per_unit_of_Capital/Labour_per_unit_of_Capital;
+ Consumption_per_unit_of_Labour=Consumption_per_unit_of_Capital/Labour_per_unit_of_Capital;
+
+ % Compute steady state of the endogenous variables.
+ L=1/(1+Consumption_per_unit_of_Labour/((1-alpha)*theta/(1-theta)*Output_per_unit_of_Labour^(1-psi)));
+ c=Consumption_per_unit_of_Labour*L;
+ k=L/Labour_per_unit_of_Capital;
+ y=Output_per_unit_of_Capital*k;
+end;
+
+steady;
+
+// Declare pre-announced shocks
+shocks;
+ var epsilon;
+ periods 5, 15;
+ values -0.1, -0.1;
+end;
+
+perfect_foresight_setup(periods=300);
+perfect_foresight_solver;
+
+// Declare unexpected shock (after first simulation!)
+oo_.exo_simul(11, 1) = 0.1; // Period 10 has index 11!
+
+saved_endo = oo_.endo_simul(:, 1:9); // Save periods 0 to 8
+saved_exo = oo_.exo_simul(1:9, :);
+oo_.endo_simul = oo_.endo_simul(:, 10:end); // Keep periods 9 to 301
+oo_.exo_simul = oo_.exo_simul(10:end, :);
+
+periods 291;
+perfect_foresight_solver;
+
+// Combine the two simulations
+oo_.endo_simul = [ saved_endo oo_.endo_simul ];
+oo_.exo_simul = [ saved_exo; oo_.exo_simul ];
+
+rplot epsilon;
+rplot c;

deterministic-examples/rbcii.mod

+var k, y, L, c, i, A, a, mu;
+varexo epsilon;
+parameters beta, theta, tau, alpha, psi, delta, rho, Astar, sigma;
+
+beta    =  0.990;
+theta   =  0.357;
+tau     =  2.000;
+alpha   =  0.450;
+psi     = -2.500;
+delta   =  0.020;
+rho     =  0.998;
+Astar   =  1.000;
+sigma   =  0.100;
+
+model;
+ a = rho*a(-1) + sigma*epsilon;
+ A = Astar*exp(a);
+ (c^theta*(1-L)^(1-theta))^(1-tau)/c - mu = beta*((c(+1)^theta*(1-L(+1))^(1-theta))^(1-tau)/c(+1)*(alpha*(y(+1)/k)^(1-psi)+1-delta)-mu(+1)*(1-delta));
+ ((1-theta)/theta)*(c/(1-L)) - (1-alpha)*(y/L)^(1-psi);
+ y = A*(alpha*(k(-1)^psi)+(1-alpha)*(L^psi))^(1/psi);
+ k = y-c+(1-delta)*k(-1);
+ i = k-(1-delta)*k(-1);
+
+[ mcp = 'i > 0' ]
+ mu = 0;
+end;
+
+steady_state_model;
+ a=0;
+ mu=0;
+ A=Astar;
+    
+ // Steady state ratios 
+ Output_per_unit_of_Capital=((1/beta-1+delta)/alpha)^(1/(1-psi));
+ Consumption_per_unit_of_Capital=Output_per_unit_of_Capital-delta;
+ Labour_per_unit_of_Capital=(((Output_per_unit_of_Capital/A)^psi-alpha)/(1-alpha))^(1/psi);
+ Output_per_unit_of_Labour=Output_per_unit_of_Capital/Labour_per_unit_of_Capital;
+ Consumption_per_unit_of_Labour=Consumption_per_unit_of_Capital/Labour_per_unit_of_Capital;
+
+ L=1/(1+Consumption_per_unit_of_Labour/((1-alpha)*theta/(1-theta)*Output_per_unit_of_Labour^(1-psi)));
+ c=Consumption_per_unit_of_Labour*L;
+ k=L/Labour_per_unit_of_Capital;
+ y=Output_per_unit_of_Capital*k;
+ i=delta*k;
+end;
+
+shocks;
+ var epsilon;
+ periods 10;
+ values -1;
+end;
+
+steady;
+
+perfect_foresight_setup(periods=400);
+perfect_foresight_solver(lmmcp, maxit=200);
+
+n = 40;
+
+figure(2);
+subplot(3,2,1); plot(1:n,A(1:n)); title('A');
+subplot(3,2,2); plot(2:n,y(2:n)); title('y');
+subplot(3,2,3); plot(2:n,L(2:n)); title('L');
+subplot(3,2,4); plot(1:n,k(1:n)); title('k');
+subplot(3,2,5); plot(2:n,c(2:n)); title('c');
+subplot(3,2,6); plot(2:n, y(2:n)-c(2:n)); title('i');
+

exercise-and-solutions/nk3.mod

+// 3-equations NK model (without ZLB): stochastic case
+
+var y     // GDP
+    pie   // Inflation
+    r     // Nominal interest rate
+    eps_a // Productivity shock (opposite of cost-push shock)
+    eps_b // Confidence shock (shock to rate of preference for the present)
+    ;
+
+varexo eta_r // Monetary policy error
+       eta_a // Innovation to productivity shock
+       eta_b // Innovation to confidence shock
+       ;
+
+parameters sigma   // Relative risk aversion
+           rho_r   // Inertia of monetary policy
+           rho_pie // Reaction of monetary policy to inflation
+           rho_y   // Reaction of monetary policy to output gap
+           beta    // Discount factor
+           rho_a   // Autocorrelation of productivity shock
+           rho_b   // Autocorrelation of confidence shock
+           kappa   // Weight of output gap in Phillips curve
+           delta   // Weight of inflation expectations in Phillips curve
+           ;
+
+sigma = 2;
+rho_r = 0.96;
+rho_pie = 1.68;
+rho_y = 0.1;
+beta = 0.99;
+rho_a = 0.82;
+rho_b = 0.85;
+delta = 0.67;
+kappa = 0.1;
+
+model;
+ // Euler equation
+ beta*r/pie(+1)*eps_b(+1)*y(+1)^(-sigma) = eps_b*y^(-sigma);
+
+ // NK Phillips curve
+ log(pie/STEADY_STATE(pie)) = delta*log(pie(+1)/STEADY_STATE(pie)) + kappa*log(y/STEADY_STATE(y)) - log(eps_a);
+
+ // Taylor rule (in multiplicative form)
+ log(r) = rho_r*log(r(-1)) + (1-rho_r)*(log(STEADY_STATE(pie)/beta) + rho_pie*log(pie(-1)/STEADY_STATE(pie)) + rho_y*log(y/STEADY_STATE(y))) + eta_r;
+
+ // Stochastic shocks (AR(1) processes)
+ log(eps_a) = rho_a*log(eps_a(-1)) + eta_a;
+ log(eps_b) = rho_b*log(eps_b(-1)) + eta_b;
+end;
+
+shocks;
+ var eta_r; stderr 0.01;
+ var eta_a; stderr 0.08;
+ var eta_b; stderr 0.05;
+end;
+
+steady_state_model;
+ eps_a = 1;
+ eps_b = 1;
+ pie = 1.005; // Steady state inflation is not determined by the model, any value will work here
+ y = 1; // Same applies to GDP
+ r = pie/beta;
+end;
+
+steady;
+check;
+stoch_simul(order=1);

exercise-and-solutions/nk3_zlb_det.mod

+// 3-equations NK model with ZLB: perfect foresight case (question 4)
+
+var y     // GDP
+    pie   // Inflation
+    r     // Nominal interest rate
+    eps_a // Productivity shock (opposite of cost-push shock)
+    eps_b // Confidence shock (shock to rate of preference for the present)
+    ;
+
+varexo eta_r // Monetary policy error
+       eta_a // Innovation to productivity shock
+       eta_b // Innovation to confidence shock
+       ;
+
+parameters sigma   // Relative risk aversion
+           rho_r   // Inertia of monetary policy
+           rho_pie // Reaction of monetary policy to inflation
+           rho_y   // Reaction of monetary policy to output gap
+           beta    // Discount factor
+           rho_a   // Autocorrelation of productivity shock
+           rho_b   // Autocorrelation of confidence shock
+           kappa   // Weight of output gap in Phillips curve
+           delta   // Weight of inflation expectations in Phillips curve
+           ;
+
+sigma = 2;
+rho_r = 0.96;
+rho_pie = 1.68;
+rho_y = 0.1;
+beta = 0.99;
+rho_a = 0.82;
+rho_b = 0.85;
+delta = 0.67;
+kappa = 0.1;
+
+model;
+ // Euler equation
+ beta*r/pie(+1)*eps_b(+1)*y(+1)^(-sigma) = eps_b*y^(-sigma);
+
+ // NK Phillips curve
+ log(pie/STEADY_STATE(pie)) = delta*log(pie(+1)/STEADY_STATE(pie)) + kappa*log(y/STEADY_STATE(y)) - log(eps_a);
+
+ // Taylor rule (in multiplicative form)
+ log(r) = max(rho_r*log(r(-1)) + (1-rho_r)*(log(STEADY_STATE(pie)/beta) + rho_pie*log(pie(-1)/STEADY_STATE(pie)) + rho_y*log(y/STEADY_STATE(y))) + eta_r, 0);
+
+ // Stochastic shocks (AR(1) processes)
+ log(eps_a) = rho_a*log(eps_a(-1)) + eta_a;
+ log(eps_b) = rho_b*log(eps_b(-1)) + eta_b;
+end;
+
+steady_state_model;
+ eps_a = 1;
+ eps_b = 1;
+ pie = 1.005; // Steady state inflation is not determined by the model, any value will work here
+ y = 1; // Same applies to GDP
+ r = pie/beta;
+end;
+
+steady;
+
+shocks;
+ var eta_a;
+ periods 1;
+ values 0.08;
+end;
+
+perfect_foresight_setup(periods=40);