var c_H, c_F, cC, o_C, ii, ppi_H, wr, o_H, y_H, mrs, pr_H, pr_O, rer, l, y_GDP, 
	X, M, y_S, O, pr_M, bstar, rr, ppi, istar, delttae, alphha_H, psi, ppistar, 
	prstar_O, cstar, vv, zetta;

varexo epsilon_istar, epsilon_O, epsilon_S, epsilon_zetta, epsilon_alphha, epsilon_ppistar, epsilon_psi, epsilon_cstar, epsilon_vv;

%Unused varexo: epsilon_delttae, epsilon_alphha_H,  epsilon_prstar_O, epsilon_delttas;

varobs y_GDP ii ppi_H rer delttae; % pr_O % wr l O;

parameters gammma thetta etta hh caprhho phhi_H betta xi_H alphha phhi_L eps_L 
	eps_H xi_L omegga sigmma_L deltta c_Hy_Hsteady cstar_Hy_Hsteady ettastar cy_steady
	NXY_steady cstar_Hx_steady mX_steady oX_steady c_Fx_steady y_Sx_steady c_Fm_steady OM_steady o_CO_steady
	o_HO_steady g_Y p_XxeB rhho_i rhho_istar rhho_ppistar omegga_ppi omegga_y omegga_delttae rhho_v rhho_O 
	rhho_S rhho_a_H rhho_alphha rhho_zetta rhho_psi rhho_ppi rhho_cstar rhho_vv;
	
gammma=.4;
thetta=1;
etta=.15;
hh=.5;
caprhho=.001;
phhi_H=.75;
betta=.99;
xi_H=.5;
alphha=.4;
phhi_L=.75;
eps_L=9;
eps_H=11;
xi_L=.5;
omegga=.1;
sigmma_L=1;
deltta=.1;
c_Hy_Hsteady=.54;
cstar_Hy_Hsteady=.36;
ettastar=.1;
cy_steady=1.2;
NXY_steady=-.25;
cstar_Hx_steady=.7;
y_Sx_steady=.3;
c_Fm_steady=.75;
OM_steady=.25;
o_CO_steady=.66;
o_HO_steady=.34;
cstar_Hx_steady=1.06;
mX_steady=1.74;
oX_steady=.44;
c_Fx_steady=1.29;
g_Y=3.5;
p_XxeB=1;
rhho_i=.75;
omegga_ppi=.75;
omegga_y=.7;
omegga_delttae=.70;
rhho_v=.7;
rhho_O=.7;
rhho_S=.7;
rhho_a_H=.7;
rhho_alphha=.7;
rhho_zetta=.7;
rhho_psi=.7;
rhho_ppi=.7;
rhho_cstar=.7;
rhho_vv=.7;
rhho_istar=.7;
rhho_ppistar=.7;

model;
#A_H=eps_H/(eps_H-1);
#v_L=betta*phhi_L;
%#zetta=((1-v_L)*(1-phhi_L))/(1+sigmma_L*eps_L);
#chii=1/((1+ppistar)*(1+g_Y));
#delttas=y_S-y_S(-1);
%Note: Forward Looking Variables
% Equation 4: cC(+1), ppi(+1)
% Equation 5: delttae(+1) -> try changing this to contemporaneous ("End of time period"), because it is the only variable that does not also have a lag in the same equation.
% Equation 7: ppi(+1)
% Equation 11: wr(+1), ppi_H(+1)

//Block 1: Aggregate Demand
c_H = (1-gammma)*(thetta-etta)*rer - (thetta*(1-gammma)+gammma*etta)*pr_H+cC; //Demand for home goods

c_F = -(thetta*gammma+etta*(1-gammma))*rer - gammma*(thetta-etta)*pr_H+cC; //Demand for foreign goods

o_C = -etta*pr_O+cC; //Oil demand

cC = (1/(1+hh))*cC(+1)+(1/(1+hh))*cC(-1)-((1-hh)/(1+hh))*(ii-ppi(+1)); //Euler Equation

ii = istar + expectation(-1)(delttae)+caprhho*bstar; 

%istar=rhho_istar*istar(-1)+epsilon_istar;

//Block 2: Aggregate Supply and Inflation
ppi_H = ((1-phhi_H)*(1-betta*phhi_H))/(phhi_H*(1+betta*xi_H))*((1-alphha)*wr+alphha*pr_O
	-alphha_H-pr_H)+betta/(1+betta*xi_H)*ppi_H(+1)+xi_H/(1+betta*xi_H)*ppi_H(-1);

o_H-l = omegga*(wr-pr_O);

y_H = alphha_H + (1-alphha)*l+alphha*o_H;

%alphha_H = rhho_a_H*a_H(-1)+epsilon_a_H;

((1+v_L*phhi_L+sigmma_L*eps_L*(phhi_L+v_L))/(1+sigmma_L*eps_L))*wr-phhi_L*wr(-1)*v_L*expectation(0)(wr(+1))
		=(((1-v_L)*(1-phhi_L))/(1+sigmma_L*eps_L))*mrs-(phhi_L+v_L*xi_L)*ppi+phhi_L*xi_L*ppi(-1)+v_L*expectation(0)(ppi(+1))+zetta;

%zetta=rhho_zetta*zetta(-1)+epsilon_zetta;

mrs=sigmma_L*l+1/(1-hh)*cC-hh/(1-hh)*cC(-1);

//Block 3: Relative Prices
pr_H=pr_H(-1)+ppi_H-ppi;

pr_O=rer+prstar_O+psi;

%prstar_O=rhho_O*prstar_O(-1)+epsilon_O;

%psi=rhho_psi*psi(-1)+epsilon_psi;

rer=rer(-1)+delttae+ppistar-ppi;

%ppistar=rhho_ppi*ppistar(-1)+epsilon_ppistar;

0=deltta*prstar_O+(1-deltta)*gammma*pr_H+(1-deltta)*(1-gammma)*rer;

//Block 4: Aggregate Equilibrium
y_H=(c_Hy_Hsteady)*c_H+(cstar_Hy_Hsteady)*cstar-ettastar
	*(cstar_Hy_Hsteady)*(pr_H-rer);

	y_GDP=(cy_steady)*cC+(NXY_steady)*(X-M);

X=-ettastar*(cstar_Hx_steady)*(pr_H-rer)+(cstar_Hx_steady)
	*cstar+(y_Sx_steady)*y_S;

	%y_S=rhho_S*y_S(-1)+epsilon_S;

%cstar=rhho_cstar*cstar(-1)+epsilon_cstar;

M=(c_Fm_steady)*c_F+(OM_steady)*O;

O=(o_CO_steady)*o_C+(o_HO_steady)*o_H;

pr_M=(c_Fm_steady)*rer+(OM_steady)*prstar_O;

(1-caprhho)*betta*bstar=betta*istar+chii*bstar(-1)+chii*X(-1)
	+(cstar_Hx_steady)*chii*pr_H+chii*(delttas-ppi)+((p_XxeB)-betta)
	*X+((p_XxeB)-betta)*pr_H-((p_XxeB)*mX_steady)*M-((p_XxeB)*c_Fx_steady)
	*rer-((p_XxeB)*oX_steady)*prstar_O;

//Block 5: Policy Rule
rr = rhho_i*rr(-1)+(1-rhho_i)*(omegga_ppi*ppi+omegga_y*y_GDP+omegga_delttae*delttae)+vv;

rr = ii-expectation(0)(ppi(+1));

//Block 6: Exogenous stochastic processes;
istar = rhho_i*istar(-1)+epsilon_istar;

alphha_H = rhho_alphha*alphha_H(-1)+epsilon_alphha;

zetta = rhho_zetta*zetta(-1)+epsilon_zetta;

%Note: Price of oil, and the associated shocks, are part of varobs so won't show up here in estimated model.
prstar_O = rhho_O*prstar_O(-1)+epsilon_O;

psi = rhho_psi*psi(-1)+epsilon_psi;

ppistar=rhho_ppistar*ppistar(-1)+epsilon_ppistar;

y_S=rhho_S*y_S(-1)+epsilon_S;

cstar=rhho_cstar*cstar(-1)+epsilon_cstar;

%Note: Monetary shocks are part of varobs so won't show up here in esimated model.
vv=rhho_vv*vv(-1)+epsilon_vv;
	
end;

estimated_params;
%Consumers
hh, normal_pdf, .5, .25;
thetta, inv_gamma_pdf, 1, .3; %inv_gamma_pdf
etta, inv_gamma_pdf, .15, .05; %inv_gamma_pdf
ettastar, normal_pdf, .1, .01;
deltta, gamma_pdf, .1, .01;
sigmma_L, normal_pdf, 1, .3;
caprhho, inv_gamma_pdf, .75, .2; %inv_gamma_pdf

%Producers
omegga, inv_gamma_pdf, .1, .5; %inv_gamma_pdf
alphha, gamma_pdf, .3, .1;

%Rigidity
phhi_H, gamma_pdf, .75, .05;
phhi_L, gamma_pdf, .75, .05;
xi_H, gamma_pdf, .1, .01;
xi_L, gamma_pdf, .1, .01;

%Monetary Policy
rhho_i, gamma_pdf, .75, .2;
omegga_ppi, normal_pdf, .75, .15;
omegga_y, normal_pdf, .75, .15; 
omegga_delttae, inv_gamma_pdf, .5, .15; %inv_gamma_pdf

%Shocks
rhho_a_H, gamma_pdf, .75, .25;
rhho_S, gamma_pdf, .75, .25;
rhho_cstar, gamma_pdf, .75, .25;
rhho_istar, gamma_pdf, .75, .25;
rhho_ppistar, gamma_pdf, .75, .25;
rhho_zetta, gamma_pdf, .75, .25;
rhho_psi, gamma_pdf, .75, .25;

end;

initval;
c_H=0;
c_F=0;
cC=0;
o_C=0;
ii=0;
ppi_H=0;
wr=0;
o_H=0;
y_H=0;
mrs=0;
rer=0;
pr_H=0;
pr_O=0;
l=0;
y_GDP=0; 
X=0;
M=0;
y_S=0; 
O=0;
pr_M=0;
bstar=0;
rr=0;

ppi=0;
istar=0;
delttae=0;
alphha_H=0;
ppistar=0;
psi=0;
prstar_O=0;
cstar=0;
%delttas=0;
vv=0;
end;

maxit_=100;
steady;

check;

model_diagnostics;
model_info;

shocks;
var epsilon_ppistar = 1;
var epsilon_istar = 1;
var epsilon_zetta = 1;
var epsilon_alphha = 1;
var epsilon_psi = 1;
%%var epsilon_O = 1;
%var epsilon_cstar = 1;
%var epsilon_S = 1;
%%var epsilon_vv = 1;
end;

estimation (datafile=varobs3, nobs=52, conf_sig=90, mh_replic=60000, mh_nblocks=15, mh_drop=.3, mh_jscale=.025, mh_init_scale=.05, mode_compute=6, mode_check, bayesian_irf, forecast=20);
%Note: Make sure the data matches the model transformation, i.e. if the model is log-linear, make sure the data is log-linearized (i.e. log-differenced).
%Note: mh_jscale needs to be tuned to get an acceptance rate between 0.2 and 0.4.  See pg. 51 of 
%the Dynare User Guide for more details.
%Note: Want to consider using mode_file, mode_compute, mode_check, load_mh_file,

identification;
dynare_sensitivity(identification=1, morris=2);