Sorry for not replying earlier. I was away for a few days. My problem came from the use of some temporary calculations involving some of the instruments in the parameters’ block. Eliminating those and following the example from M. Monfils solved my issue. I have a unique optimal policy for this benchmark model. I have implemented an identical approach in other models with success so thanks for your help and patience.
I am however puzzled with a new thing. Csolve returns a steadty state capital calculation that I know is wrong for one of my models. I checked it by running dynare with instruments in the parameters block and running my steady state code with fsolve. Note that the csolve discrepancy is only true for low values of some of the instruments. Any idea of what may cause this? I am thinking this is may be down to the algorithm used by Csolve. I have included the code below.
Many thanks,
Script:
A = 1;
theta1=0.3;
theta2=0.6;
theta3=0.1;
deltak=0.1;
deltag=0.1;
gamma=3;
b=1.5;
beta=0.99;
deltaka=0;
alpha=0.05;
eta = 1/(gamma - theta2);
m = (theta2/(b*gamma))^eta;
rho1 = theta1*( 1 + theta2*eta)-1;
rho3 = theta3*( 1 + theta2*eta );
rhotaul = theta2*eta;
rhotfp = 1 + theta2*eta;
omega = (1/(theta1*m^theta2))^(1/rho3);
zeta1 = theta1*eta - rho1/rho3*theta3*eta;
zeta2 = eta - rhotfp/rho3*eta*theta3;
zetataul = eta*(1-theta3*rhotaul/rho3);
tauc=0.02;
taud=0.9;
taul=0.02;
tauk=0.02;
TFP = A;
ka = (alpha/((1-tauc)*(1/beta-1+deltaka)))^(1/(1-alpha))
%Fsolve call
x0=[2];
options = optimset('Algorithm','levenberg-marquardt','MaxFunEvals',2000,'Tolx',10e-12,'Tolfun',10e-12);
[SS, fval, exitflag,ouput]=fsolve(@(x) divd4fsolve(x,beta,tauk,tauc,taud,taul,theta1,theta2,theta3,b,gamma,deltag,deltak,eta,m,omega,rho1,rho3,rhotfp,rhotaul,zeta1,zeta2,zetataul,TFP,alpha,deltaka,ka),x0,options);
k=SS(1)
%Csolve call
k = find_kdivd4(2.3,beta,taud,tauc,taul,tauk,theta1,theta2,theta3,b,gamma,deltag,deltak,eta,m,omega,rho1,rho3,rhotfp,rhotaul,zeta1,zeta2,zetataul,TFP,alpha,deltaka,ka);
Fsolve:
function g = divd4fsolve(x,beta,tauk,tauc,taud,taul,theta1,theta2,theta3,b,gamma,deltag,deltak,eta,m,omega,rho1,rho3,rhotfp,rhotaul,zeta1,zeta2,zetataul,TFP,alpha,deltaka,ka);
k= x(1);
g(1) = deltag*omega*((1/beta-1)/(1-tauk)+deltak)^(1/rho3)*(k^(-rho1/rho3)*TFP^(-rhotfp/rho3))*(1-taul)^(-rhotaul/rho3)...
-( tauk*theta1 + theta2*taul + theta3*(taud + tauc*(1-taud)))*TFP*k^theta1*...
(m*(1-taul)^zetataul*((1/beta-1)/(1-tauk)+deltak)^(theta3*eta/rho3)*omega^(theta3*eta)*k^zeta1*TFP^zeta2)^theta2*...
(omega*((1/beta-1)/(1-tauk)+deltak)^(1/rho3)*(k^(-rho1/rho3)*TFP^(-rhotfp/rho3))*(1-taul)^(-rhotaul/rho3))^theta3...
+ tauk*deltak*k...
- taud*alpha*ka^alpha...
+ (taud + tauc*(1-taud))*deltaka*ka;
end
Csolve:
function k = find_kdivd4(k0,beta,tauk,tauc,taud,taul,theta1,theta2,theta3,b,gamma,deltag,deltak,eta,m,omega,rho1,rho3,rhotfp,rhotaul,zeta1,zeta2,zetataul,TFP,alpha,deltaka,ka)
k = csolve(@divd_ss4,k0,],1e-15,100,beta,tauk,tauc,taud,taul,theta1,theta2,theta3,b,gamma,deltag,deltak,eta,m,omega,rho1,rho3,rhotfp,rhotaul,zeta1,zeta2,zetataul,TFP,alpha,deltaka,ka);
function g = divd_ss4(k,beta,tauk,tauc,taud,taul,theta1,theta2,theta3,b,gamma,deltag,deltak,eta,m,omega,rho1,rho3,rhotfp,rhotaul,zeta1,zeta2,zetataul,TFP,alpha,deltaka,ka);
g = deltag*omega*((1/beta-1)/(1-tauk)+deltak)^(1/rho3)*(k^(-rho1/rho3)*TFP^(-rhotfp/rho3))*(1-taul)^(-rhotaul/rho3)...
-( tauk*theta1 + theta2*taul + theta3*(taud + tauc*(1-taud)))*TFP*k^theta1*...
(m*(1-taul)^zetataul*((1/beta-1)/(1-tauk)+deltak)^(theta3*eta/rho3)*omega^(theta3*eta)*k^zeta1*TFP^zeta2)^theta2*...
(omega*((1/beta-1)/(1-tauk)+deltak)^(1/rho3)*(k^(-rho1/rho3)*TFP^(-rhotfp/rho3))*(1-taul)^(-rhotaul/rho3))^theta3...
+ tauk*deltak*k...
- taud*alpha*ka^alpha...
+ (taud + tauc*(1-taud))*deltaka*ka;
end
end