Bytecode does not enforce positivity constraint on irreversible investment model
In the following RBC model with irreversible investment, the positivity constraint is not enforced if simulated with block+bytecode. On the contrary, if block+bytecode is removed, the classic deterministic solver gives the right result.
var k, y, l, c, i, A, a, expterm, 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(differentiate_forward_vars,block,bytecode,cutoff=0);
a = rho*a(-1) + sigma*epsilon;
A = Astar*exp(a);
mu = max(0,(((c^theta)*((1-l)^(1-theta)))^(1-tau))/c - expterm(1)+beta*mu(1)*(1-delta));
(i<=0)*(k - (1-delta)*k(-1)) + (i>0)*((((c^theta)*((1-l)^(1-theta)))^(1-tau))/c - expterm(1)+beta*mu(1)*(1-delta)) = 0;
expterm = beta*((((c^theta)*((1-l)^(1-theta)))^(1-tau))/c)*(alpha*((y/k(-1))^(1-psi))+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);
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;
expterm=beta*(c^theta*(1-l)^(1-theta))^(1-tau)/c*(1+alpha*(y/k)^(1-psi)-delta);
end;
shocks;
var epsilon;
periods 10;
values -1;
end;
steady;
options_.maxit_ = 200;
simul(periods=400);
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');