Fix treatment of complex elements in perfect foresight simulations
The mod-file
%----------------------------------------------------------------
% define variables
%----------------------------------------------------------------
var y c k n r w invest G lambda q phi; %I excluded ghat
varexo eps_g;
%----------------------------------------------------------------
% define parameters
%----------------------------------------------------------------
parameters beta eta b mu delta alpha gamma tau_k tau_n kappa gshare;
%----------------------------------------------------------------
% set parameter values
%----------------------------------------------------------------
beta=1.03^(-1/4);
b=0;
mu=0;
delta=0.021;
alpha = 0.34;
gamma=1.005;
tau_k=0;
tau_n=0;
kappa=0;
gshare= 0.2038;
%----------------------------------------------------------------
% enter model equations
%----------------------------------------------------------------
model;
1/(c-c(-1)*b/gamma)-beta*b*1/(c(+1)*gamma-b*c)=lambda;
q=beta*lambda(+1)/(lambda*gamma)*((1-delta)*q(+1)+(1-tau_k)*r+delta*tau_k);
1=q*(1-kappa/2*(invest/invest(-1)*gamma-gamma)^2-kappa*(invest/invest(-1)*gamma-gamma)*invest*gamma/invest(-1))+beta*q(+1)*lambda(+1)/(lambda*gamma)*kappa*(invest(+1)/invest*gamma-gamma)*(invest(+1)*gamma/invest)^2;
lambda*w*(1-tau_n)=eta/((1-n)^mu);
c+invest=(1-tau_n)*w*n+(1-tau_k)*r*k(-1)+delta*tau_k*k(-1)-phi;
G=tau_n*w*n+tau_k*(r-delta)*k(-1)+phi;
k(-1)=w/r*n*alpha/(1-alpha);
k*gamma=(1-delta)*k(-1)+(1-kappa/2*(invest/invest(-1)*gamma-gamma)^2)*invest;
y=k(-1)^alpha*n^(1-alpha);
y=c+invest+G;
G=gshare*steady_state(y)*exp(eps_g);
end;
%----------------------------------------------------------------
% set steady state values and calibrate the model to steady state labor of 0.24, i.e. compute the corresponding steady state values
% and the labor disutility parameter by hand
%---------------------------------------------------------------
steady_state_model;
n=0.24;
q=1;
r=(q*gamma/beta-(1-delta)*q-delta*tau_k)/(1-tau_k);
k=n*(1/alpha*r)^(1/(alpha-1));
y=r*k/alpha;
G=gshare*y;
c=r*k/alpha-(gamma+delta-1)*k-G;
w=(c+(gamma+delta-1-r)*k+G)/n;
invest=(gamma+delta-1)*k;
phi=G-tau_n*w*n-tau_k*(r-delta)*k;
lambda=(1/c)*((1-b/gamma)^(-1)-beta*b*(gamma-b)^(-1));
eta=lambda*w*(1-tau_n)*(1-n)^(mu);
end;
shocks;
var eps_g;
periods 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50;
values -0.008075285 0.006460228 0.024225855 0.046836652 0.062987222 0.074292621 0.077522735 0.080752848 0.082367905 0.080752848 0.080752848 0.079137791 0.077522735 0.075907678 0.071062507 0.067832393 0.064602279 0.061372165 0.058142051 0.054911937 0.051681823 0.050066766 0.046836652 0.043606538 0.041991481 0.040376424 0.038761367 0.03714631 0.035531253 0.032301139 0.030686082 0.029071025 0.027455968 0.027455968 0.025840912 0.025840912 0.024225855 0.024225855 0.022610798 0.022610798 0.020995741 0.020995741 0.019380684 0.019380684 0.017765627 0.017765627 0.01615057 0.01615057 0.014535513 0.012920456;
end;
%----------------------------------------------------------------
% check the starting values for the steady state
%---------------------------------------------------------------
resid(1);
%----------------------------------------------------------------
% compute steady state given the starting values
%---------------------------------------------------------------
steady;
%----------------------------------------------------------------
% check Blanchard-Kahn-conditions
%---------------------------------------------------------------
check;
simul(periods=100,maxit=1000,stack_solve_algo=6);
terminates with imaginary parts in the results. Judging from the size of imaginary parts (maximum in the range of 1e-13) this is due to numerical inaccuracies.
With stack_solve_algo=6
, the message is Perfect foresight solution found
while with stack_solve_algo=0
, it fails with
Simulation terminated with imaginary parts in the residuals or endogenous variables.
We should aim at making this consistent. My proposal would be to discard any imaginary parts smaller than 1e-10 in perfect foresight simulations and return the cleaned values as a solution.
Sidenote: in the attached mod-file, capital explodes despite the model being Blanchard-Kahn stable. This is something regularly reported by users in the forum in the presence of endogenous stocks. IS there a good explanation for this?