%endogeneous variables
var c h pnyoc  pie mc piem mcm lamda f wnew w piec y pnyom pnyo rh blbar rl pl q bbarld bs r x m cdo z d g rs bl blstar
    hstar piecstar mcstar lamdastar fstar wnewstar wstar ystar rhstar blbarstar rlstar plstar bsstar rstar zstar dstar gstar cstar
    dy_obs dw_obs dh_obs pie_obs r_obs rl_obs dystar_obs dwstar_obs dhstar_obs piecstar_obs rstar_obs rlstar_obs q_obs;

%exogeneous variables
varexo ebs ebl ebbarl ebbarlstar ebsstar eblstar
       omegaz omegablbar omegabs omegad omegag omegazstar omegablbarstar omegadstar omegabsstar omegagstar
       uy uw uh upie ur url uq uystar uwstar uhstar upiestar urstar urlstar omegar omegarstar;

%//parameters
parameters eta theta phir phipiec phiy rhoz rhoblbar rhod rhobs rhog
           beta ky ksiy km ksim b sigmac ksiw kappaw phi xnyo
           kystar ksiystar ksiwstar phistar xnyostar thetastar phirstar phicstar phiystar
           rhozstar rhoblbarstar rhodstar rhobsstar rhogstar
           lamdam hss lamdax blyss gstarss bsstarystarss xyss gss
           kappawstar sigmacstar bstar
           cstaryss blstarystarss lamdaw lamdawstar
           zss blystarss lamday lamdaystar hstarss zstarss bsss ebsss ebbarlss ebsstarss ebbarlstarss;

beta = 0.9;            %lit_rev
zss=1;                 %lit_rev
hss=1;                 %lit_rev
gstarss = 0.18;        %lit_rev
hstarss=1.00;          %lit_rev
zstarss=1.00;          %lit_rev
lamdaystar = 1.20;     %lit_rev
lamdawstar =1.05;      %lit_rev
gss = 0.17;            %lit_rev
lamdam = 1.20;         %lit_rev
lamdax = 1.20;         %lit_rev
lamday = 1.20;         %lit_rev
lamdaw = 1.05;         %lit_rev
cstaryss =0.1;         %tekito
%bsyss =0.3;
bsstarystarss =0.1;    %tekito
xyss =0.2;             %tekito
bsss=0.1;              %tekito

blyss = 0.2;
theta = 0.75;
xnyo = 10;
eta = 1.5;
phipiec = 1.5;
phiy = 0.13;
phir = 0.75;
sigmac = 1.5;
phi = 2;
b = 0.7;
ky = 0.5;
ksiy = 0.5;
km = 0.5;
ksim = 0.5;
kappaw = 0.5;
ksiw = 0.5;
rhoz = 0.75;
rhoblbar = 0.75;
rhod = 0.75;
rhobs = 0.75;
rhog = 0.75;

blstarystarss = 0.2;
thetastar = 0.75;
xnyostar = 10;
phicstar = 1.5;
phiystar = 0.13;
phirstar = 0.75;
sigmacstar = 1.5;
phistar = 2;
bstar = 0.7;
kystar = 0.5;
ksiystar = 0.5;
kappawstar = 0.5;
ksiwstar = 0.5;
rhozstar = 0.75;
rhoblbarstar = 0.75;
rhodstar = 0.75;
rhobsstar = 0.75;
rhogstar = 0.75;

ebsss=1;       
ebbarlss=1;
ebsstarss=1;
ebbarlstarss=1;
%blstaryss=0.1;

model(linear);
#kappa = 1;       
#kappastar = 1;
#piecss=1;
#piess=1;
#piemss=1;

#rhss=1/beta;
#rsss=1/(ebsss*beta);
#rlss=(rsss*ebsss)/ebbarlss;
#rsstarss=1/(ebsstarss*beta);
#rlstarss=(rsstarss*ebsstarss)/ebbarlstarss;

#yss=zss*hss;
#mcmss=1/lamdam;
#pnyomss=lamdam;
#mcxss=1/lamdax;
#pnyoxss=1;
#pnyoss=pnyoxss/lamdax;
%#T=lamdam*lamdax;

#yblbarss=(beta/(1-beta*kappa))*blyss;
#blss=((1-beta*kappa)/beta)*yblbarss;
#ystarblbarstarss=(beta/(1-beta*kappastar))*blstarystarss;
#blstarss=((1-beta*kappastar)/beta)*(ystarblbarstarss);
#blbarss=blss/(rlss-1);

#blstaryss=blstarss/yss;
#blbarstaryss=blbarss/yss;


#teta=1-(blss/blstarss);
#tera=teta/(1-teta);

#bbarlstar=(beta/(1-beta*kappastar))*(1+tera)*blss;
#delta=bsss/(2*blss);
#bsyss=(yblbarss+blbarstaryss)*(1/delta);
#deltastar = (bsstarystarss/blstarystarss)*(1/2);

#alpha = (((1/(lamdam^(1-eta)))-((1/(lamdam*lamdax))^(1-eta))))/(1-((1/(lamdam*lamdax))^(1-eta)));
#cdoyss=((pnyoss*eta)/(1-alpha))*(1-(xyss)-gss);
#myss=xyss-(1-rlstarss)*(bbarlstar/yss);
#wss=(pnyoss*zss)*lamday;
#css=(1-gss)*zss*hss;
#lamdass=(1-beta*b)/(((1-b)^sigmac)*(css^sigmac));
#psi=((1/lamdaw)*lamdass*wss)/(hss^phi);

#ystarss=zstarss*hstarss;
#piecstarss=1;
#mcstarss=1/lamdaystar;
#wstarss=zstarss/lamdaystar;
#cstarss=(1-gstarss)*zstarss*hstarss;
#lamdastarss=1/((cstarss)^sigmac)*(((1-beta*b)*d)/((1-b)^sigmac));

%#eta_mc=lamdam/(lamdam-1);
%#gamma_cmc=((1-alpha)*((eta_mc-1)/eta_mc)^(1-eta)+alpha)^(1/(1-eta));
%#gamma_mcc=1/gamma_cmc;



%DOMESTIC EQUATIONS
%//inflation(1)
pie = (beta/(1+beta*ky))*pie(+1) + (ky/(1+beta*ky))*pie(-1)+ ((1-ksiy)*(1-beta*ksiy)/(ksiy*(1+beta*ky)))*mc;

%//marginal cost(2)
mc = w - pnyo - z;

%//import price inflation(3)
piem = (beta/(1+beta*km))*piem(+1) + (km/(1+beta*km))*piem(-1) + ((1-ksim)*(1-beta*ksim)/(ksim*(1+beta*km)))*mcm;

%//importers mc(4)
mcm = q - pnyom;

%//marginal utility of comsumption(5)
lamda = (1/(1-b*beta)) * ((d-b*beta*d(+1)) - ((sigmac/(1-b)) * ((1+beta*(b^2))*c-b*beta*c(+1)-b*c(-1))));

%//comsumption euler equation(6)
lamda = lamda(+1) + (lamda(+2)-lamda(+1))+ rh - piec(+1);

%//wage fillips curve1(7)
f = (1-beta*ksiw) * (((-1/(lamdaw-1))*wnew)+lamda+((lamdaw/(lamdaw-1))*w)+h) + (beta*ksiw)* ((1/(lamdaw-1))*(piec(+1)-(kappaw*piec)+(wnew(+2)-wnew(+1)))+f(+1));

%//wage fc2(8)
f = (1-beta*ksiw) * (d+(1+phi)*(lamdaw/(lamdaw-1))*(w-wnew)+(1+phi)*h) + (beta*ksiw) * ((1+phi)*(lamdaw/(lamdaw-1))*(piec(+1)-(kappaw*piec)+(wnew(+2)-wnew(+1)))+f(+1));

%//wage fc3(9)
wnew-w = (ksiw/(1-ksiw)) * (piec-(kappaw*piec(-1))+w-w(-1));

%//cpi inflation(10)
%piec-pnyoc(-1) = (1-alpha)*pie + alpha*pnyomss*(1-eta)*(piem+pnyom(-1)-pnyo(-1));
pnyoc = (1-alpha)*pie(+1) + alpha*pnyomss*(1-eta)*(piem(+1)+pnyom-pnyo) - piec(+1);
%pnyoc=alpha*(gamma_mcc)^(1-eta)*pnyom;

%//production function(11)
y = h + z;

%//goods market clearing (12)
y = cdoyss*cdo + xyss*x + gss*(q+y);

%//import relative prices (13)
pnyom = pnyom(-1) + piem - piec;

%//domestic relative prices (14)
pnyo = pnyo(-1) + pie - piec;

%//household effective rate (15)
rh = rs + ebs + xnyo * delta * (((blss/(bl+((1-teta)*blstarss)))*bl)+((((1-teta)*blstarss)/(blss+((1-teta)*blstarss)))*blstar)-bs);

%//long-term debt value (16)
blbar = pl + bl;

%//long-term interest rate (17)
rl(+1)+pl(+1)-pl = rh - ebl + (xnyo*(rlss - kappa))*(((blss/(blss+(1 - teta)*blstarss))*bl)+((((1 - teta)*blstarss)/(blss+(1 - teta)*blstarss))*blstar)-bs);
%rl+pl-pl(-1) = rh(-1) - ebl(-1)+ (xnyo*(rlss - kappa))*((blss/(blss+(1 - teta)*blstarss)*bl(-1))+((((1 - teta)*blstarss)/(blss+(1 - teta)*blstarss))*blstar(-1))-bs(-1));

%//UIP condition (18)
%rlstar(+1)-piecstar(+1)+plstar(+1)-plstar+q(+1)-q = rh - eblstar + ((xnyo*(rlstarss-kappastar))/rlstarss) * ((blss/(blss+(1-teta)*blstarss)*bl + (((1-teta)*blstarss)/(blss+(1-teta)*blstarss))*blstar - bs));
rlstar-piecstar+plstar-plstar(-1)-q-q(-1) = rh(-1) - eblstar(-1) + ((xnyo*(rlstarss-kappastar))/rlstarss) * ((blss/(blss+(1-teta)*blstarss)*bl(-1) + (((1-teta)*blstarss)/(blss+(1-teta)*blstarss))*blstar(-1) - bs(-1)));

%//foreign long-term debt value (19)
bbarld = q + plstar + blstar;

%//long-term asset price (20)
pl = -(1/(1-beta*kappa))*rl;

%//short-term debt (21)
bsyss*(bs-ebs) + blyss*(blbar-ebbarl) + bsyss*theta*bs(-1) = (1/beta) * bsyss * (bs(-1)+r(-1)-piec) + (1/beta) * blyss * (blbar(-1)+rl+(pl-pl(-1))-piec) + pnyoss*gss*(pnyo+g+y);

%//taylor rule (22)
r = phir*r(-1) + (1-phir)*((phipiec*piec)+phiy*y) + omegar;

%//exports (23)
xyss*x = myss * (q+m) + blstaryss * (bbarld-ebbarlstar) - (1/beta)*blstaryss*(bbarld(-1)+rlstar-piecstar+(plstar-plstar(-1))+(q-q(-1)));

%//imports (24)
m = -(eta*pnyom)+c;

%//domestic consumption (25)
cdo = -(eta*pnyo)+c;

%//TFP (26)
z = (rhoz*z(-1)) + omegaz ;

%//long-term debt risk premium (27)
blbar = rhoblbar * blbar(-1) + omegablbar;

%//discount factor (28)
d = rhod * d(-1) + omegad;

%//short-term debt risk premium (29)
bs = rhobs * bs(-1) + omegabs;

%//government spending (30)
g = rhog * g(-1) + omegag;


%%//FOREIGN EQUATIONS

%//Inflation (31)
piecstar = (beta/(1+beta*kystar)) * piecstar(+1) + (kystar/(1+beta*kystar)) * piecstar(-1) + (((1-ksiystar)*(1-beta*ksiystar))/(ksiystar*(1+beta*kystar))) * mcstar;

%//marginal cost (32)
mcstar = wstar - zstar;

%//marginal utility of consumption
lamdastar = (1/(1-bstar*beta)) * ((dstar-bstar*beta*dstar(+1)) - (sigmacstar/(1-bstar)) * ((1+beta*(bstar^2))*cstar - bstar*beta*cstar(+1)-bstar*cstar(-1)));

%//mcu (33)
lamdastar = lamdastar(+1) + (dstar(+2) - dstar(+1)) + rhstar - piecstar(+1);

%//wage fc1 (34)
fstar = (1-beta*ksiwstar) * ((-1/(lamdawstar-1)) * wnewstar + lamdastar + (lamdawstar/(lamdawstar-1)) * wstar + hstar) + (beta * ksiwstar) * ((1/(lamdawstar-1)) * (piecstar(+1) - (kappawstar * piecstar) + (wnewstar(+2)-wnewstar(+1))) + fstar(+1));

%//wage fc2 (35)
fstar = (1-beta*ksiwstar) * (dstar + (1+phistar) * (lamdawstar/(lamdawstar-1)) * (wstar - wnewstar) + (1+phistar) * hstar) + (beta * ksiwstar) * ((1+phistar) * (lamdawstar/(lamdawstar-1)) * (piecstar(+1) - (kappawstar * piecstar) + (wnewstar(+2)-wnewstar(+1))) + fstar(+1));

%//wage fc3 (36)
wnewstar - wstar = (ksiwstar/(1-ksiwstar)) * (piecstar - (kappawstar * piecstar(-1)) + wstar - wstar(-1));

%//production function (37)
ystar = hstar + zstar;

%//goods market clearing (38)
ystar = cstaryss * cstar + gstarss * (gstar + ystar);

%//household effective rate (39)
rhstar = rstar + ebsstar + xnyostar * deltastar * (blstar - bsstar);

%//long-term debt value (40)
blbarstar = plstar + blstar;

%//long-term interst rate (41)
%rlstar(+1) + (plstar(+2)-plstar(+1)) = rhstar - ebbarlstar + xnyostar * (rlstarss - kappastar) * (blstar - bsstar);
rlstar + plstar - plstar(-1) = rhstar(-1) - ebbarlstar(-1) + xnyostar * (rlstarss - kappastar) * (blstar(-1) - bsstar(-1));

%//long-term asset price (42)
plstar = (-1/(1-beta*kappastar)) * rlstar;

%//short-term debt (43)
bsstarystarss* (bsstar - ebsstar) + blstarystarss * (blbarstar - ebbarlstar) + bsstarystarss * thetastar * bsstar(-1) = gstarss * (gstar + ystar) * (1/beta) * bsstarystarss * (bsstar(-1) + rstar(-1) - piecstar) + (1/beta) * blstarystarss * (blstar(-1) + rlstar + (plstar(+1)-plstar) - piecstar);

%//taylor rule (44)
rstar = phirstar * rstar(-1) + (1-phirstar) * (phicstar*piecstar+phiystar*ystar) + omegarstar;

%//TFP (45)
zstar = rhozstar * zstar(-1) + omegazstar;

%//long-term debt risk premium (46)
blbarstar = rhoblbarstar * blbarstar(-1) +  omegablbarstar;

%//dicount factor (47)
dstar = rhodstar * dstar(-1) + omegadstar;

%//short-term debt risk premium (48)
bsstar = rhobsstar * bsstar(-1) + omegabsstar;

%//governmetn spending (49)
gstar = rhogstar * gstar(-1) + omegagstar;

%//observation EQUATIONS
dy_obs = y + uy;
dw_obs = w + uw;
dh_obs = h + uh;
pie_obs = pie + upie;
r_obs = r + ur;
rl_obs = rl + url;
dystar_obs = ystar + uystar;
dwstar_obs = wstar + uwstar;
dhstar_obs = hstar + uhstar;
piecstar_obs = piecstar + upiestar;
rstar_obs = rstar + urstar;
rlstar_obs = rlstar + urlstar;
q_obs = q + uq;

end;


NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = deblank(M_.endo_names(ii,:));
  eval([varname '=0;']);
end


%check(qz_zero_threshold=1e-10);
model_diagnostics;

estimated_params;
%blbsblstyss, normal_pdf, 0.4, 0.1;
blyss, normal_pdf, 0.2, 0.01;
theta, normal_pdf, 0.75, 0.25;
xnyo, normal_pdf, 10, 1;
eta, normal_pdf, 1.5, 0.1;
phipiec, normal_pdf, 1.5, 0.1;
phiy, normal_pdf, 0.13, 0.05;
phir, beta_pdf, 0.75, 0.05;
sigmac, normal_pdf, 1.5, 0.1;
phi, normal_pdf, 2, 0.1;
b, beta_pdf, 0.7, 0.1;
ky, beta_pdf, 0.5, 0.15;
ksiy, beta_pdf, 0.5, 0.1;
km, beta_pdf, 0.5, 0.15;
ksim, beta_pdf, 0.5, 0.05;
kappaw, beta_pdf, 0.5, 0.15;
ksiw, beta_pdf, 0.5, 0.1;
rhoz, beta_pdf, 0.75, 0.1;
rhoblbar, beta_pdf, 0.75, 0.1;
rhod, beta_pdf, 0.75, 0.1;
rhobs, beta_pdf, 0.75, 0.1;
rhog, beta_pdf, 0.75, 0.1;
stderr uy, inv_gamma_pdf, 0.1, inf;
stderr uw, inv_gamma_pdf, 0.1, inf;
stderr uh, inv_gamma_pdf, 0.1, inf;
stderr upie, inv_gamma_pdf, 0.1, inf;
stderr ur, inv_gamma_pdf, 0.1, inf;
stderr url, inv_gamma_pdf, 0.1, inf;
stderr uq, inv_gamma_pdf, 0.1, inf;

blstarystarss, normal_pdf, 0.2, 0.01;
thetastar, normal_pdf, 0.75, 0.25;
xnyostar, normal_pdf, 10, 1;
phicstar, normal_pdf, 1.5, 0.1;
phiystar, normal_pdf, 0.13, 0.05;
phirstar, beta_pdf, 0.75, 0.1;
sigmacstar, normal_pdf, 1.5, 0.1;
phistar, normal_pdf, 2, 0.1;
bstar, beta_pdf, 0.7, 0.1;
kystar, beta_pdf, 0.5, 0.15;
ksiystar, beta_pdf, 0.5, 0.1;
kappawstar, beta_pdf, 0.5, 0.15;
ksiwstar, beta_pdf, 0.5, 0.1;
rhozstar, beta_pdf, 0.75, 0.1;
rhoblbarstar, beta_pdf, 0.75, 0.1;
rhodstar, beta_pdf, 0.75, 0.1;
rhobsstar, beta_pdf, 0.75, 0.1;
rhogstar, beta_pdf, 0.75, 0.1;
stderr uystar, inv_gamma_pdf, 0.1, inf;
stderr uwstar, inv_gamma_pdf, 0.1, inf;
stderr uhstar, inv_gamma_pdf, 0.1, inf;
stderr upiestar, inv_gamma_pdf, 0.1, inf;
stderr urstar, inv_gamma_pdf, 0.1, inf;
stderr urlstar, inv_gamma_pdf, 0.1, inf;

end;

varobs dy_obs dw_obs dh_obs pie_obs r_obs rl_obs dystar_obs dwstar_obs dhstar_obs piecstar_obs rstar_obs rlstar_obs q_obs;

estimated_params_init(use_calibration);
end;

estimation(datafile = jpdat, plot_priors=0, mode_check, mh_replic=500000, mh_nblocks=2, mh_drop=0.5, mh_jscale=0.5, bayesian_irf);

write_latex_dynamic_model;