From f78ac9f5ac5563fdc30abadd6e77388bbdb5fd08 Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx.de>
Date: Wed, 26 May 2021 12:35:27 +0200
Subject: [PATCH] EP with MCP: add test case

---
 tests/Makefile.am      |  1 +
 tests/ep/rbcii_MCP.mod | 98 ++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 99 insertions(+)
 create mode 100644 tests/ep/rbcii_MCP.mod

diff --git a/tests/Makefile.am b/tests/Makefile.am
index e76ce7fc1..52052c2e2 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -324,6 +324,7 @@ MODFILES = \
 	ep/linearmodel0.mod \
 	ep/linearmodel1.mod \
 	ep/rbc_bytecode.mod \
+	ep/rbcii_MCP.mod \
 	stochastic_simulations/example1_noprint.mod \
 	stochastic-backward-models/solow_cd.mod \
 	stochastic-backward-models/solow_ces.mod \
diff --git a/tests/ep/rbcii_MCP.mod b/tests/ep/rbcii_MCP.mod
new file mode 100644
index 000000000..698896d2d
--- /dev/null
+++ b/tests/ep/rbcii_MCP.mod
@@ -0,0 +1,98 @@
+% RBC model  with  irreversible  investment  constraint, implemented using
+% MCP tag
+
+var k, y, L, c, A, a, mu, i;
+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 -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);
+i=y-c;
+[ mcp = 'i > 0' ]
+mu = 0;
+
+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;
+i=delta*k;
+mu=0;
+end;
+
+steady;
+
+%% Perfect  foresight
+% negative  technology shock drives the  economy  to  the  bound
+%shocks;
+%var epsilon;
+%periods 1:20;
+%values -0.5;
+%end;
+%perfect_foresight_setup(periods=300);
+%perfect_foresight_solver(lmmcp, maxit=200);
+%figure
+%subplot(2,2,1)
+%plot(oo_.endo_simul(8,:)); axis  tight
+%title('investment')
+%subplot(2,2,2)
+%plot(oo_.endo_simul(2,:)); axis  tight
+%title('output')
+%subplot(2,2,3)
+%plot(oo_.endo_simul(3,:)); axis  tight
+%title('hours')
+%subplot(2,2,4)
+%plot(oo_.endo_simul(4,:)); axis  tight
+%title('consumption')
+
+
+shocks;
+var epsilon;
+stderr 0.10;
+end;
+
+extended_path(periods=300,lmmcp);
+
+if any(oo_.endo_simul(strmatch('i',M_.endo_names,'exact'),:)<-1e-6)
+    error('lmmcp tag did not work')
+end
+% // Plot 100 first periods of  the simulations
+% figure
+% subplot(2,2,1)
+% plot(oo_.endo_simul(strmatch('i',M_.endo_names,'exact'),1:101)); axis  tight
+% title('investment')
+% subplot(2,2,2)
+% plot(oo_.endo_simul(strmatch('y',M_.endo_names,'exact'),1:101)); axis  tight
+% title('output')
+% subplot(2,2,3)
+% plot(oo_.endo_simul(strmatch('L',M_.endo_names,'exact'),1:101)); axis  tight
+% title('hours')
+% subplot(2,2,4)
+% plot(oo_.endo_simul(strmatch('c',M_.endo_names,'exact'),1:101)); axis  tight
+% title('consumption')