close all;

var 
c l1 c1 cg1 ct1 cn1 ch1 cf1 B1 R w   
FI FI1 ER EB1 rp
pt  
I pi It In Ih If 
Pr_n  yn kn  l_n xn gn in  pkn  Psi_n Psi1_n rn pn  
Pr_t  yt kt  l_t xt gt it  pkt  Psi_t Psi1_t rt ph  
Pr_x  yx kx  l_x xx gx ix  pkx  Psi_x Psi1_x rx px
Pr_g  yg kg  l_g xg gg ig  pkg  Psi_g Psi1_g rg pg
e y_eu pf x_eu g_eu
GDP
;
predetermined_variables kn kt kx kg EB1 B1;
varexo e_rp;

parameters alpha gamma alpha_g gamma_g theta a_t a_n a_ch a_it a_ih An At Ax Ag phi beta delta psi0 psi1 psi2 g_b ;
    alpha   = 0.3;
    gamma   = 0.4;
    alpha_g   = 0.4;
    gamma_g   = 0.3;
    theta   = 0.25;
    a_it    = 0.38;
    a_ih    =0.31;
    a_t     =0.23;
    a_n     =0.73;
    a_ch    =0.68;
    At       = 2.3;
    An       = 1.7;
    Ax      =2.1;
    Ag      =3;
    phi     = 0.9;
    beta    = 0.99;
    delta   = 0.025;
    psi2    =  -0.005/delta+1;
    psi0    =  (psi2-1)/psi2*delta;
    psi1    =  delta^(1-psi2)/psi2;   
    g_b=1;



model; 

%%%%%%%Non-tradable%%%%%%%
c=c1;
log(yn)=log(An)+alpha*log(kn)+gamma*log(l_n)+theta*log(gn)+(1-alpha-gamma-theta)*log(xn);
kn(+1)=(1-delta+Psi_n)*kn;
gamma*pn*yn=w*l_n;
(1-alpha-gamma-theta)*yn*pn=px*xn;
theta*yn*pn=pg*gn;
pkn*Psi1_n=pi;
(c(+1)*pkn)/(c*pkn(+1))=beta*(1+rn(+1)-delta+Psi_n-in(+1)*Psi1_n/kn(+1));
kn*pkn*rn=alpha*yn*pn;
Psi_n=psi0+psi1*(in/kn)^psi2;
Psi1_n=psi1*psi2*(in/kn)^(psi2-1);
Pr_n=pn*yn-pi*in-w*l_n-px*xn-pg*gn;

%%%%%%%Tradable%%%%%%%
log(yt)=log(At)+alpha*log(kt)+gamma*log(l_t)+theta*log(gt)+(1-alpha-gamma-theta)*log(xt);
kt(+1)=(1-delta+Psi_t)*kt;
gamma*ph*yt=w*l_t;
(1-alpha-gamma-theta)*yt*ph=px*xt;
theta*yt*ph=pg*gt;
pkt*Psi1_t=pi;
(c(+1)*pkt)/(c*pkt(+1))=beta*(1+rt(+1)-delta+Psi_t-it(+1)*Psi1_t/kt(+1));
kt*pkt*rt=alpha*yt*ph;
Psi_t=psi0+psi1*(it/kt)^psi2;
Psi1_t=psi1*psi2*(it/kt)^(psi2-1);
Pr_t=ph*yt-pi*it-w*l_t-px*xt-pg*gt;

%%%%%%%Oil%%%%%%%
log(yx)=log(Ax)+alpha*log(kx)+gamma*log(l_x)+theta*log(gx)+(1-alpha-gamma-theta)*log(xx);
kx(+1)=(1-delta+Psi_x)*kx;
gamma*px*yx=w*l_x;
(1-alpha-gamma-theta)*yx=xx;
theta*yx*px=pg*gx;
pkx*Psi1_x=pi;
(c(+1)*pkx)/(c*pkx(+1))=beta*(1+rx(+1)-delta+Psi_x-ix(+1)*Psi1_x/kx(+1));
kx*pkx*rx=alpha*yx*px;
Psi_x=psi0+psi1*(ix/kx)^psi2;
Psi1_x=psi1*psi2*(ix/kx)^(psi2-1);
Pr_x=px*yx-pi*ix-w*l_x-px*xx-pg*gx;


%%%%%%%Gas%%%%%%%
log(yg)=log(Ag)+alpha_g*log(kg)+gamma_g*log(l_g)+theta*log(gg)+(1-alpha_g-gamma_g-theta)*log(xg);
kg(+1)=(1-delta+Psi_g)*kg;
gamma_g*pg*yg=w*l_g;
(1-alpha_g-gamma_g-theta)*yg*pg=px*xg;
theta*yg=gg;
pkg*Psi1_g=pi;
(c(+1)*pkg)/(c*pkg(+1))=beta*(1+rg(+1)-delta+Psi_g-ig(+1)*Psi1_g/kg(+1));
kg*pkg*rg=alpha_g*yg*pg;
Psi_g=psi0+psi1*(ig/kg)^psi2;
Psi1_g=psi1*psi2*(ig/kg)^(psi2-1);
Pr_g=pg*yg-pi*ig-w*l_g-px*xg-pg*gg;

%%%%%%%Consumer1%%%%%%%
c1+B1(+1)/R+e*(EB1(+1)/((1-FI)*ER))=B1+e*EB1+w*l1+(Pr_x+Pr_n+Pr_g+Pr_t);
ln(c1)=a_n*ln(cn1)+a_t*ln(ct1)+(1-a_n-a_t)*ln(cg1);
ln(ct1)=a_ch*ln(ch1)+(1-a_ch)*ln(cf1);
pn/pt=(a_n*ct1)/(a_t*cn1);
pg/pt=((1-a_n-a_t)*ct1)/(a_t*cg1);
ph/(e*pf)=(a_ch*cf1)/((1-a_ch)*ch1);
ln(pt)=a_ch*ln(ph/a_ch)+(1-a_ch)*ln(e*pf/(1-a_ch));

phi*c1=w*(1-l1);
c1(+1)/c1=beta*R;
1=beta*(ER*e(+1)*c1/(c1(+1)*e)*(((1-FI)^2)/(1-FI+EB1(+1)*FI1)));
FI=g_b*(exp((e*EB1(+1))/GDP-STEADY_STATE(EB1))-1)-rp;
FI1=e*g_b*exp((e*EB1(+1))/GDP-STEADY_STATE(EB1))/GDP;



%%%%%%%Investment%%%%%
ln(I)=a_it*ln(It)+(1-a_it)*ln(In);
ln(It)=a_ih*ln(Ih)+(1-a_ih)*ln(If);
e*pf*If=a_it*(1-a_ih)*pi*I;
ph*Ih=a_it*a_ih*pi*I; 
pn*In=(1-a_it)*pi*I; 

%%%%%%%Market-clearing%%%%%%%
yn=cn1+In;
yt=ch1+Ih+y_eu; 

I=in+ix+ig+it;
yx=xx+xn+xg+xt+x_eu;
yg=cg1+gx+gn+gg+g_eu;
%pg=0.6;
%px=0.6;
l1=l_n+l_x+l_g+l_t;
B1=0;


pg*g_eu+px*x_eu+ph*y_eu-e*pf*(cf1+If)=e*(EB1(+1)/ER-EB1);

GDP=pn*yn+ph*yt+px*yx+pg*yg;

%%%%%%%Exogeneous%%%%%%%
y_eu=0.1;
pf=0.7;
g_eu=0.09;
x_eu=0.12;
rp=0.6*rp(-1)+e_rp;

ER=1;

end;

initval;

c      	=	0.504661	;
l1     	=	0.384442	;
c1     	=	0.504661	;
cg1    	=	0.10445	;
ct1    	=	0.165511	;
cn1    	=	0.781686	;
ch1    	=	0.226582	;
cf1    	=	0.0849137	;
B1     	=	0	;
R      	=	1.0101	;
w      	=	0.737859	;
FI     	=	-0.362503	;
FI1    	=	0.397497	;
ER     	=	1	;
EB1    	=	1.19585	;
rp     	=	0	;
pt     	=	0.529657	;
I      	=	0.211765	;
pi     	=	0.827936	;
It     	=	0.116526	;
In     	=	0.305393	;
Ih     	=	0.0785039	;
If     	=	0.139152	;
Pr_n   	=	0.0334049	;
yn     	=	1.08708	;
kn     	=	3.99438	;
l_n    	=	0.209764	;
xn     	=	0.0671431	;
gn     	=	0.662733	;
in     	=	0.0998594	;
pkn    	=	0.827936	;
Psi_n  	=	0.025	;
Psi1_n 	=	1	;
rn     	=	0.035101	;
pn     	=	0.355946	;
Pr_t   	=	0.00920064	;
yt     	=	0.405086	;
kt     	=	1.10016	;
l_t    	=	0.0577747	;
xt     	=	0.018493	;
gt     	=	0.182535	;
it     	=	0.027504	;
pkt    	=	0.827936	;
Psi_t  	=	0.025	;
Psi1_t 	=	1	;
rt     	=	0.035101	;
ph     	=	0.26309	;
Pr_x   	=	0.00625232	;
yx     	=	0.25134	;
kx     	=	0.747617	;
l_x    	=	0.039261	;
xx     	=	0.012567	;
gx     	=	0.124042	;
ix     	=	0.0186904	;
pkx    	=	0.827936	;
Psi_x  	=	0.025	;
Psi1_x 	=	1	;
rx     	=	0.035101	;
px     	=	0.288146	;
Pr_g   	=	0.0219815	;
yg     	=	1.3083	;
kg     	=	2.62843	;
l_g    	=	0.0776427	;
xg     	=	0.0331367	;
gg     	=	0.327075	;
ig     	=	0.0657108	;
pkg    	=	0.827936	;
Psi_g  	=	0.025	;
Psi1_g 	=	1	;
rg     	=	0.035101	;
pg     	=	0.145964	;
e      	=	0.471949	;
y_eu   	=	0.1	;
pf     	=	0.7	;
x_eu   	=	0.12	;
g_eu   	=	0.09	;
GDP    	=	0.756902	;


end;


shocks;
   var e_rp = 1;
end;
steady(maxit = 1000, solve_algo = 0);
check;
model_diagnostics;
%stoch_simul(order = 2, irf=15);
