// ***Preamble section***
// Define the endogenous variables
var k a l c y w r g;
// Define the exogenous variables (shocks)
varexo epsg epsa;
// Define the parameters of the model
parameters alpha rhoa rhog rho lstar gamma delta gystar;
// Set the values of the parameters to match Romer’s calibration
alpha = 1/3;
rhoa  = .95;
rhog  = .95;
rho = 0.01;
lstar = 1/3;
gamma = 0.005;
delta = 0.025;
gystar = 0.2;

// ***Model section***
model;

// Equation (19) from above. Expected future value of x is x(+1)
-c = r(+1) - c(+1);

// Equation (17) from above. Lags of x denoted by x(-1)
(1+gamma) * k = (rho + gamma + delta)/alpha  * y
	+ (1-delta) * k(-1)
	+ (gamma + delta - (rho + gamma + delta)*(1-gystar)/alpha) * c
	- gystar*(rho+gamma+delta)/alpha * g;

// Equation (13) from above.
y = alpha * k(-1) + (1-alpha) * (a + l);

// YOU WILL NEED TO ADD THREE EQUATIONS HERE

// Equations (7) and (8) from above
a = rhoa * a(-1) + epsa;
g = rhog * g(-1) + epsg;

end;
// *** Steady-state section *** In the steady state all deviations are zero
initval;
c = 0;
l = 0;
k = 0;
y = 0;
g = 0;
a = 0;
w = 0;
r = 0;
end;

steady;
check;

// *** Shocks section ***
shocks;
var epsa;
periods 1;
values .01;
end;

// *** Computation section ***
simul(periods=200);

// *** Some MATLAB plot commands similar to those below 
subplot(3,3,1);
plot(k);
title('K');
axis([0 40 -0.002 .01]);
subplot(3,3,2);
plot(a);
title('A');
axis([0 40 -0.002 .01]);
subplot(3,3,3);
plot(l);
title('L');
axis([0 40 -0.002 .01]);
subplot(3,3,4);
plot(c);
title('C');
axis([0 40 -0.002 .01]);
subplot(3,3,5);
plot(y);
title('Y');
axis([0 40 -0.002 .01]);
subplot(3,3,6);
plot(w);
title('W');
axis([0 40 -0.002 .01]);
subplot(3,3,7);
plot(r);
title('R');
axis([0 40 -0.002 .01]);
subplot(3,3,8);
plot(g);
title('G');
axis([0 40 -0.002 .01]);