var K I N M B i_f Y P D G T C w h qk theta V e q p GAMMAn GAMMAu 
Pstar Pf Df SC BS lambda; %28

%varexo ;

parameters delta s gi rho phik chi psi alpha A etal sigman sigmau 
csi eta csistar etastar omegabar epsilon phib beta Ystar ti istar tauc tauw tauf ;

delta = 0.006;
s = 0.024;
gi = 0.25;
rho = 0.56;
phik = 7;
chi = 0.941;
psi = 0.6;
alpha = 0.32;
A = 1;
etal = 1/0.56;
sigman = 18.31;
sigmau = 10.96;
csi = 0.8;
eta = 1.01;
csistar = 0.993;
etastar = 1.01;
omegabar = 0.13;
epsilon = 0.6;
phib = 0.0019;
beta = 0.99;
Ystar = 0.91;
ti = 0.09;
istar = 0.02;

tauc = 0.22;
tauw = 0.13;
tauf = 0.34;


options = optimset('Display','off','TolFun',1e-10);
x0 = [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0];
x = fsolve(@steadystate, x0, options, csistar, etastar, Ystar, csi, eta, omegabar, chi, psi, beta, s, alpha, tauf, sigman, sigmau, etal, tauw, tauc, epsilon, A, gi, rho, delta, ti)

Pfss = x(1);
Yss = x(2);
Pss = x(3);
thetass = x(4); 
Nss = x(5); 
wss = x(6);
hss = x(7);
Css = x(8);
ess = x(9); 
pss = x(10); 
qss = x(11); 
Iss = x(12); 
Vss = x(13);
k = x(14);
Kss = x(15);






model;
% 1
K = (1 - delta)*K(-1) + I;
% 2
N = (1 - s)*N(-1) + M;
% 3
B = (1 + i_f)*B(-1) + Y/P - D;
% 4  
P*G = gi*Y;
% 5
P*G + P*T = tauc*P*C + (tauf + tauw)*N(-1)*w*h - (1 - tauw)*rho*w*h*(1 - N(-1));
% 6
P*T = ti*Y;
% 7
qk = 1 + phik*(I - delta*K(-1))/K(-1);
% 8
theta = V/(e*(1 - N(-1)));
% 9
q = chi*theta^(psi - 1);
% 10
p = M/(e*(1 - N(-1)));
% 11
Y = A*K(-1)^alpha*(N(-1)*h)^(1 - alpha);
% 12
M = chi*V^psi*(e*(1 - N(-1)))^(1 - psi);
% 13
GAMMAn = -sigman*h^(1 + etal)/(1 + etal);
% 14
GAMMAu = -sigmau*e^(1 + etal)/(1 + etal);
% 15
Y = csi*P^eta*D + (1 - csistar)*Pstar^(etastar)*Ystar;
% 16
D = C + I + omegabar*V + phik/2*(I - delta*K(-1))^2/K(-1) + G;
% 17
P = (csi + (1 - csi)*Pf^(1 - eta))^(1/(1 - eta));
% 18
Pstar = (csistar*Pf^(1 - etastar) + (1 - csistar))^(1/(1 - etastar));
% 19
Df = (1 - csi)*(Pf/P)^(-eta)*D;
% 20
SC = omegabar*((1 - s)/q*(1 - (1 + tauf)/(1 + tauf(+1))*(1 - tauw(+1))/(1 - tauw)) + e*theta*((1 + tauf)/(1 + tauf(+1))*(1 - tauw(+1))/(1 - tauw)));
% 21
BS = 1/P*(1 - alpha)*Y/N(-1) + SC;
% 22
%1/C = (1 + tauc)*P*lambda;
% 23
sigman*h^(1 + etal) = (1 - tauw)/((1 + tauc)*(1 + tauf))*(1 - alpha)*Y/N(-1)*1/(P*C);
% 24
sigmau*e^(1 + etal) = (1 - tauw)/((1 + tauc)*(1 + tauf))*(1 - epsilon)/epsilon*omegabar*e*theta*1/C;
% 25
w*h = (1 - epsilon)/(1 + tauf)*P*BS + epsilon*(rho*w*h + (GAMMAu - GAMMAn)/((1 - tauw)*lambda));
% 26
1 + i_f = 1 + istar - phib*P*B(-1)/Y;
% 27   
omegabar/q = beta*(P(+1)*lambda(+1)/(P*lambda)*(1/P(+1)*(1 - alpha)*Y(+1)/N + (1 - s)*omegabar/q(+1) - (1 + tauf(+1))*w(+1)*h(+1)));
% 28
qk = beta*(P(+1)*lambda(+1)/(P*lambda)*(1/P(+1)*alpha*Y(+1)/K + qk(+1) - delta + phik/2*((I(+1) - delta*K)/K)^2));
% 29
1 = beta*(P(+1)*lambda(+1)/(P*lambda)*(1 + istar(+1) - phib*P(+1)*B(+1)/Y(+1)));
%
end;

initval;

%tauc = 0.22;
%tauw = 0.13;
%tauf = 0.34;

K = Kss;
I = Iss; 
N = Nss;

M = 0.39*0.27*1.001; 

B = 0; 

i_f = istar; 

Y = Yss; 
P = Pss; 

D = 0.86/1.91; 
G = gi*0.86/1.91; 
T = ti*0.86/1.91; 

C = Css; 
w = wss; 
h = hss; 

qk = 1; 

theta = thetass; 
V = Vss; 
e = ess; 
q = qss;
p = pss; 

GAMMAn = -sigman*0.28^(1 + etal)/(1 + etal); 
GAMMAu = -sigmau*0.27^(1 + etal)/(1 + etal);

Pstar = 23; 

Pf = Pfss; 

Df = (1 - csi)*(27.07/1.91)^(-eta)*0.86/1.91; 
SC = omegabar*((1 - s)/1.64*(1 - (1 + tauf)/(1 + tauf)*(1 - tauw)/(1 - tauw)) + 0.27*7.21*((1 + tauf)/(1 + tauf)*(1 - tauw)/(1 - tauw)));
BS = 1/1.91*(1 - alpha)*0.86/1.001 + SC;
lambda = 1/(0.28*1.91)*1/(1 + tauc); 

end;

steady;

 

simul(periods=100, stack_solve_algo=6, maxit=100);