%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 0. Housekeeping
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. Variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

var  c, c1, c2, y, y1, h, h1, q, Bbar, Y, w, l, l1, L, b, b1; 
varexo theta;
parameters beta, gamma, v, u, rho, phi, eta, sigma, R, A; 

beta=0.95;
gamma=0.90;
u=0.36;
v=0.10;
phi=1.05;
rho=0.95;
eta=3;
sigma=0.014;
A=1;
R=1/gamma;



model;
  c1 + q*h1 + R*b1(-1)+w*l1 = y1 + q*h1(-1) + b1;
  c + q*h + R*b(-1) +w*l= y+ q*h(-1) + b;
  c+c1+c2+(R*(b1(-1)+b(-1)))+(q*(h1+h))= Y+(q*(h1(-1)+h(-1)))+(b1+b);
  q(-1)= (beta*q) + (beta*u*((y1/h1(-1))));
  q(-1)= (gamma*v*((y/h(-1)))) + (gamma*q) + ((beta-gamma)*theta*v*((y/h(-1))));
  R*b(-1)=theta*y;
  h+h1 = 1;
  b+b1=0;
  l+l1=L;
  L=(1/w)^(1+1/eta);
  1/c2=(1/gamma)*((L^(1/eta))/w);
  c2=w*L;
  y = A*(h(-1)^v)*(l(-1)^(1-v));
  y1 = Bbar*(h1(-1)^u)^(l1(-1)^(1-u));
  Y=y+y1;
  phi=A/Bbar;

end;

 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3. Computation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
initval;
 h = 0.276;
 h1 = 1-h;
 L=(1/gamma)^(1/((-(eta^2)-eta-1)/((eta+1)*eta)-1/(eta+1)));
 w=L^(-eta/(eta+1));
 c2=w*L;
 l=((((((1/phi)*((1-u)/(1-v))*(((h1^u)/(h^v))))^(1/v))^1/v)+1)^-1)*L;
 l1=L-l;
 q = (gamma/(1-((theta*beta) + ((1-theta)*gamma))))*A;
 Bbar=0.95;
 y = A*(h^v)*(l^(1-v));
 y1 = Bbar*(h1^u)*(l1^(1-u));
 b = theta*beta*y;
 Y=y+y1; 
 c = y + ((1-(1/beta))*b); 
 c1 = y1 + ((1-(1-1/beta))*b1);
 theta=0.5; 

end;
steady;
check;

shocks;
var theta; periods 1;
values 0.20;
end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4. Simulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

simul(periods=2100);