Hi, did you get the code to work? My steady state seemed wrong but I think I checked everything and I don’t know what to do now. Here is my code, any help will be greatly appreciated.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Can news about the future drive the business cycle?
%
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calibration and Preliminary Computation
sigma0=1;
theta0=1.4; % labor elasticity
beta0=0.985; % subjective discount factor
alpha0=0.64; % labor share in production
gamma0=0.001;
delta0 =0.025;
psi0=2.3; % only needs to be greater than 0
phi0 = 1.3; % Adjustment cost param
para2=0.025/(1/0.15-0.5);
para1= (1/0.15-1)*0.025/(1/0.15-0.5); %Util Param.
rho_a=0.95;
rho_z=0.95;
%% Definition of the Rational Expectations Model
% state variables: k0 a0 z0 x0 epsa0 epsz0
% z0 state or choice??
% choice variables: y0 c0 i0 n0 w0 r0 u0 delta_u0 lambda0 eta0 mu0
syms k0 a0 X0 epsa0 epsz0
syms k1 a1 X1 epsa1 epsz1
syms z0 c0 i0 n0 gdp0 u0 lambda0 eta0 mu0 if0
syms z1 c1 i1 n1 gdp1 u1 lambda1 eta1 mu1 if1
%arguments of the F function
x0 = [k0 a0 X0 epsa0 epsz0];
x1 = [k1 a1 X1 epsa1 epsz1];
y0 = [z0 c0 i0 n0 gdp0 u0 lambda0 eta0 mu0 if0];
y1 = [z1 c1 i1 n1 gdp1 u1 lambda1 eta1 mu1 if1];
nx=length(x0); ny=length(y0);n0=nx+ny;
% equilibrium equations
f1=(c0-(psi0*(n0^theta0)*X1))^(-sigma0)+(mu0*gamma0*(X0/c0)^(1-gamma0))-lambda0; %7
f2=(c0-(psi0*(n0^theta0)*X1))^(-sigma0)*psi0*(n0^theta0) + mu0 - beta0*mu1*(1-gamma0)*((c1/X1)^gamma0); %8
f3=(c0-(psi0*(n0^theta0)*X1))^(-sigma0)*theta0*psi0*(n0^(theta0-1))*X1 - lambda0*alpha0*a0*(u0*k0/n0)^(1-alpha0); %9
f4=a0*(u0*k0)^(1-alpha0)*n0^alpha0-c0 - i0/z0; %5
% quadratic form
% f(u0)=para1*u0+para2/2*u0^2
% phi(i/i(-1))= phi0/2*(if0/i0-1)^2
f5=lambda0*(1-alpha0)*a0*(n0/u0)^alpha0*k0^(1-alpha0)-eta0*k0*(para1+para2*u0); %10
f6=eta0-beta0*lambda1*(1-alpha0)*a1*(n1/k1)^alpha0*u1^(1-alpha0) - eta1*(1-(para1*u0+para2/2*u0^2)); %11
f7=lambda0/z0 - eta0*(1-(phi0/2*(if0/i0-1)^2)-phi0*(if0/i0-1)*(if0/i0))-beta0*eta1*phi0*(if1/if0-1)*(if1/if0)^2; %12
f8=X1-(c0^gamma0)*X0^(1-gamma0); %2
f9= k1- i0*(1-phi0/2*(if0/i0-1)^2)-(1-(para1*u0+para2/2*u0^2))*k0 ; %6
f10=gdp0-i0/z0; %4
f11=if0-i1;
% exogenous processes
f12=log(a1)-rho_a*log(a0)-log(epsa0);
f13=log(z1)-rho_z*log(z0)-log(epsz0);
f14=log(epsa1)-0.0*log(epsa0);
f15=log(epsz1)-0.0*log(epsz0);
% f function
f=[f1;f2;f3;f4;f6;f6;f7;f8;f9;f10;f11;f12;f13;f14;f15];
z= [x0 x1 y0 y1];
fl = subs(f, z, exp(z));
%% Steady State
as=1;zs=1;us=1;epsa=1;epsz=1;
nkratio=((1/beta0-1+delta0)/(1-alpha0))^(1/alpha0); %11& %12
ckratio=(nkratio^alpha0-delta0); %4
cnratio=ckratio/(nkratio); %
% Para10=(alpha0*nkratio^(alpha0-1));
% Para20=(theta0*cnratio-Para10*gamma0/((beta0*(1-gamma0)-1)))*psi0;
% N_ss=(Para10/Para20)^(1/theta0);
paran1 = theta0*psi0*cnratio/(alpha0*nkratio^(alpha0-1));
paran2 = psi0*gamma0/(beta0*(1-gamma0)-1);
N_ss=(1/(paran1-paran2))^(1/theta0);
%%#N_ss=(Para3/psi)^(1/theta);
%psi=alpha*nkratio^(alpha-1)/((theta*cnratio+((gamma/(beta*(1-gamma)-1))*alpha*nkratio^(alpha-1)))*N_ss^theta);
C_ss=N_ss*cnratio; %
K_ss=N_ss/nkratio; %
I_ss=K_ss*delta0; %
Y_ss=K_ss^(1-alpha0)*N_ss^alpha0; %
X_ss=C_ss; %
Para1=(C_ss-psi0*N_ss^theta0*X_ss)^(-sigma0); %
% Para2=(C_ss-psi0*N_ss^theta0*X_ss);
Mu_ss=Para1*psi0*N_ss^theta0/(beta0*(1-gamma0)-1); %
Lam_ss=Para1+gamma0*Mu_ss; %
% etas = Lam_ss*(1-alpha0)*Y_ss/(K_ss*(para1+para2*us));
etas =Lam_ss;
%log transformation of the steady state values
n0=log(N_ss);n1=n0;
c0=log(C_ss);c1=c0;
k0=log(K_ss);k1=k0;
i0=log(I_ss);i1=i0;
if0=log(I_ss);if1=if0;
gdp0=log(Y_ss);gdp1=gdp0;
X0=log(X_ss);X1=X0;
mu0 =log(Mu_ss);mu1=mu0;
lambda0 =log(Lam_ss);lambda1=lambda0;
eta0 =log(etas);eta1=eta0;
u0=log(us);u1=u0;
a0=log(as);a1=a0;
z0=log(zs);z1=z0;
epsa0=log(epsa);epsa1=epsa0;
epsz0=log(epsz);epsz1=epsz0;
%test the steady state
fs=double(eval(fl(:)))
if max(abs(fs))>=1e-10; error('Incorrect Approximation Point at 0'); end
%% Approximation procedure
%Compute analytical derivatives of f
[fx,fxp,fy,fyp,fypyp,fypy,fypxp,fypx,fyyp,fyy,fyxp,fyx,fxpyp,fxpy,fxpxp,fxpx,fxyp,fxy,fxxp,fxx]=anal_deriv(fl,x0,y0,x1,y1);
%Order of approximation desired
approx = 1;
%Obtain numerical derivatives of f evaluated at the steady state
num_eval
%First-order approximation
[gx,hx,exitflag] = gx_hx(nfy,nfx,nfyp,nfxp);
%% IRF
T=100;X0=1:T;
% Technology shock
x0=[0;0;0;0;1;0];
[IR]=ir(gx,hx,x0,T);
fig01=figure('Color','w','Position',[100 100 600 400]);
ind={'1' '2' '5' '7'};
name={'c' 'n' 'gdp' 'a'};
ni=length(name);
for i0=1:ni
j=eval(ind{i0});
subplot(2,2,i0), plot(X0,IR(:,j),'LineWidth',2,'MarkerSize',5);
title(name{i0},'FontSize',10);
set(gca,'xtick',0:10:T,'FontSize',10)
axis([X0(1) X0(end) ylim]);
% ticks_format('%2.0f', '%2.2f');
end
% Government spending shock
x0=[0;0;0;0;0;1];
[IR]=ir(gx,hx,x0,T);
fig02=figure('Color','w','Position',[100 100 600 400]);
X0=1:T;
ind={'1' '2' '5' '8'};
name={'c' 'n' 'gdp' 'g'};
ni=length(name);
for i0=1:ni
j=eval(ind{i0});
subplot(2,2,i0), plot(X0,IR(:,j),'LineWidth',2,'MarkerSize',5);
title(name{i0},'FontSize',10);
set(gca,'xtick',0:10:T,'FontSize',10)
axis([X0(1) X0(end) ylim]);
% ticks_format('%2.0f', '%2.2f');
end
%% Second Moments
% zx=[eye(nx);gx]; %all variables
% eta=[0 1]';
% varshock=csigea^2*eta*eta';
% [sigz,sigx]=mom(zx,hx,varshock);
%%
%end rbc.m