%% Monacelli's model with habit formation and output gap
%% Revised on 02/13/2007
%% By Daniel Beltran
%% NOTE: Run this using Dynare

%----------------------------------------------------------------
% 0. Housekeeping (close all graphic windows)
%----------------------------------------------------------------

close all;

%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------

%% Declare Endogenous Vars
var pi_h mc y z x c pi_f psi q pie int y_star int_star pi_star y_gap s mc_star y_gap_star z_star x_flex y_flex c_flex y_star_flex;

%% Declare Exogenous Vars
varexo eph epf eq ei ei_star ep_star ez ez_star;

%% Declare Parameters
parameters sigma h phi eta theta_h theta_f rho_z rho_z_star rho_i psi_pi psi_y_tilda rho_i_star psi_pi_star psi_y_tilda_star bet gam sigma_ph sigma_pf sigma_q sigma_i sigma_i_star sigma_p_star sigma_z sigma_z_star kappa_h kappa_f BB CC DD HH JJ;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------

%% Initialize parameters:
sigma=1.0;
h=0.5;
phi=1.0;
eta=1.0;
theta_h=0.75;
theta_f=0.75;
rho_z=0.96;
rho_z_star=0.96;
rho_i=0.5;
psi_pi=1.5;
psi_y_tilda=0.5;
rho_i_star=0.6;
psi_pi_star=2;
psi_y_tilda_star=0.5;
sigma_i=0.01; 
sigma_z=0.01; 
sigma_i_star=0.01; 
sigma_z_star=0.01; 
sigma_ph=0.01; 
sigma_pf=0.01; 
sigma_q=0.01; 
sigma_p_star=0.01; 

%% other parameters:
bet=0.99;
gam=0.3;
kappa_h=(1-theta_h)*(1-bet*theta_h)/theta_h;
kappa_f=(1-theta_f)*(1-bet*theta_f)/theta_f;
BB=gam*eta*(2-gam);
CC=sigma/(sigma+phi*(1-h));
DD=gam*(2-gam)*eta-(1-h)/sigma;
HH=(1-h)/sigma;
JJ=(1-h)*(1+phi)/((1-h)*phi+sigma);

%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------

%% Model Declaration
model(linear);
pi_h=kappa_h*mc+bet*pi_h(+1)+eph;
mc=phi*y-(1+phi)*z+gam*x+(sigma/(1-h))*(c-h*c(-1));
pi_f=kappa_f*psi+bet*pi_f(+1)+epf;
q=psi+(1-gam)*x;
x-x(-1)=pi_f-pi_h;
c-h*c(-1)=y_star-h*y_star(-1)+(1/sigma)*(1-h)*((1-gam)*x+psi);
int-pie(+1)-(int_star-pi_star(+1))=q(+1)-q+eq;
(1-gam)*c=y-gam*eta*(2-gam)*x-gam*eta*psi-gam*y_star;
z=rho_z*z(-1)+ez;
y_star-h*y_star(-1)=y_star(+1)-h*y_star-(1/sigma)*(1-h)*(int_star-pi_star(+1));
pi_star=bet*pi_star(+1)+kappa_h*mc_star+ep_star;
//int=rho_i*int(-1)+(1-rho_i)*(psi_pi*pie+psi_y_tilda*y_gap)+ei;
int_star=rho_i_star*int_star(-1)+(1-rho_i_star)*(psi_pi_star*pi_star+psi_y_tilda_star*y_gap_star)+ei_star;
mc_star=phi*y_star-(1+phi)*z_star+(sigma/(1-h))*(y_star-h*y_star(-1));
z_star=rho_z_star*z_star(-1)+ez_star;
psi-psi(-1)=s-s(-1)+pi_star-pi_f;
pie=pi_h+gam*(x-x(-1));
phi*y_flex-(1+phi)*z+gam*x_flex+(sigma/(1-h))*(c_flex-h*c_flex(-1));
c_flex-h*c_flex(-1)=y_star_flex-h*y_star_flex(-1)+(1/sigma)*(1-h)*((1-gam)*x_flex);
(1-gam)*c_flex=y_flex-gam*eta*(2-gam)*x_flex-gam*y_star_flex;
y_gap=y-y_flex;
phi*y_star_flex-(1+phi)*z_star+(sigma/(1-h))*(y_star_flex-h*y_star_flex(-1));
y_gap_star=y_star-y_star_flex;
end;

%----------------------------------------------------------------
% 4. Computation
%----------------------------------------------------------------

initval;
pi_h=0;
mc=0;
y=0;
z=0;
x=0;
c=0;
pi_f=0;
psi=0;
q=0;
pie=0;
int=0;
y_star=0;
int_star=0;
pi_star=0;
y_gap=0;
s=0;
mc_star=0;
y_gap_star=0;
z_star=0;
x_flex=0;
y_flex=0;
c_flex=0;
y_star_flex=0;
end;

shocks;
var eph;stderr sigma_ph; 
var epf;stderr sigma_pf; 
var eq;stderr sigma_q; 
var ei;stderr sigma_i; 
var ei_star;stderr sigma_i_star; 
var ep_star;stderr sigma_p_star; 
var ez;stderr sigma_z; 
var ez_star;stderr sigma_z_star; 
end;

steady;
%stoch_simul;

olr_inst int;
optim_weights;
pie 2; y_gap 1;
end;
olr;