var c l_s l_n d
    c_e i_e b_e k_e y_e A_e lamd_e llamd_e x mpl_e mpk_e 
    ratio xx i_s b_s k_s y_s ps A_s lamd_s1 lamd_s2  mpl_s mpk_s
    c_b lamd_b llamd_b
    is ie 
    z1 z2 inflp infl vp y sdf 
    R_d R_e R_s w_s w_n q_e q_s 
    mp 
    lny ;

varexo u_z u_zs u_ze u_mp;

parameters 
    sigma beta_h  zeta_n alpha_e alpha
    m_b beta_b beta_e m_e delta beta_s  m_s alpha_s 
    thete epsi control_xx fy_e fy_s Rs rho_r 
    fy_pi fy_y
    cs l_ss l_ns ds
    c_es i_es b_es k_es y_es A_es lamd_es llamd_es xs mpl_es mpk_es 
    ratios  i_ss b_ss k_ss y_ss pss A_ss lamd_s1s lamd_s2s  mpl_ss mpk_ss
    c_bs lamd_bs llamd_bs
    iss ies 
    z1s z2s    sdfs 
    mpi_es, mpi_ss, i_costss, i_costes, inflps, infls, vps, ys, q_es, q_ss
    R_ds R_es R_ss w_ss w_ns ;
load data2 ;
set_param_value('cs',cs);
set_param_value('l_ss',l_ss);
set_param_value('l_ns',l_ns);
set_param_value('ds',ds);
set_param_value('c_es',c_es);
set_param_value('i_es',i_es);
set_param_value('b_es',b_es);
set_param_value('k_es',k_es);
set_param_value('y_es',y_es);
set_param_value('A_es',A_es);
set_param_value('lamd_es',lamd_es);
set_param_value('llamd_es',llamd_es);
set_param_value('mpl_es',mpl_es);
set_param_value('mpk_es',mpk_es);
set_param_value('xs',xs);
set_param_value('ratios',ratios);
set_param_value('i_ss',i_ss);
set_param_value('b_ss',b_ss);
set_param_value('k_ss',k_ss);
set_param_value('y_ss',y_ss);
set_param_value('pss',pss);
set_param_value('A_ss',A_ss);
set_param_value('lamd_s1s',lamd_s1s);
set_param_value('lamd_s2s',lamd_s2s);
set_param_value('mpl_ss',mpl_ss);
set_param_value('mpk_ss',mpk_ss);
set_param_value('c_bs',c_bs);
set_param_value('lamd_bs',lamd_bs);
set_param_value('llamd_bs',llamd_bs);
set_param_value('iss',iss);
set_param_value('ies',ies);
set_param_value('z1s',z1s);
set_param_value('z2s',z2s);
set_param_value('R_ds',R_ds);
set_param_value('R_es',R_es);
set_param_value('R_ss',R_ss);
set_param_value('w_ss',w_ss);
set_param_value('w_ns',w_ns);
set_param_value('mpi_es',mpi_es);
set_param_value('mpi_ss',mpi_ss);
set_param_value('i_costss',i_costss);
set_param_value('i_costes',i_costes);
set_param_value('inflps',inflps);
set_param_value('infls',infls);
set_param_value('vps',vps);
set_param_value('ys',ys); 
set_param_value('q_ss',q_ss);
set_param_value('q_es',q_es);
set_param_value('sigma',sigma) ;
set_param_value('beta_h',beta_h) ;
set_param_value('beta_b',beta_b) ;
set_param_value('beta_e',beta_e) ;
set_param_value('beta_s',beta_s) ;
set_param_value('zeta_n',zeta_n) ;
set_param_value('alpha_e',alpha_e) ;
set_param_value('alpha',alpha) ; 
set_param_value('alpha_s',alpha_s) ;
set_param_value('m_e',m_e) ;
set_param_value('m_b',m_b) ;
set_param_value('m_s',m_s) ;
set_param_value('delta',delta) ;
set_param_value('thete',thete) ;
set_param_value('sdfs',sdfs) ;
epsi = 10 ;
control_xx = 20;
fy_e = 0.5 ;
fy_s = 0.5;
Rs = 1beta_h ;
rho_r = 0.6 ; 
fy_pi = 1.5 ;
fy_y = 0.1;




model;
%%%% household
1 =  sdf(1) R_d   infl(1) ; 
c^(-sigma)  w_n = l_n^(zeta_n) ;
c^(-sigma)  w_s = l_s^(zeta_n) ;

 
%%%%  
c_e + q_e  i_e + b_e(-1)  R_e(-1)infl = x  y_e + b_e - w_n  l_n - ps  y_s  ;
y_e = A_e  (k_e(-1)^(alpha_e)  l_n^(1 - alpha_e))^(alpha)  y_s^(1 - alpha) ;
i_e = k_e - (1 - delta)  k_e(-1) ;
b_e = m_e  q_e(1)  k_e  infl(1)  R_e ;
lamd_e = c_e ^(-sigma) ;
mpk_e = alpha_e  alpha  y_ek_e(-1) ;
mpl_e = (1- alpha_e)  alpha  y_e l_n ;
q_e lamd_e = beta_e  lamd_e(1)  ((1 - delta)  q_e(1) + x(1) mpk_e(1)) + llamd_e  m_e q_e(1)  infl(1)R_e ;
w_n  = x  mpl_e ;
ps = x  (1 - alpha)  y_e y_s ;
lamd_e - lamd_e(1)  beta_e  R_einfl(1) - llamd_e = 0 ;
 
%%%% 
ratio = exp(xx) (1 + exp(xx)) ;
xx= control_xx  log(yy(-1)) ;
y_s = A_s  k_s(-1) ^(alpha_s)  l_s^(1 - alpha_s) ;
i_s = k_s - (1 - delta)  k_s(-1) ;
b_s = m_s  q_s(1)  k_s  infl(1)  R_s ;
mpk_s =  alpha_s  y_s  k_s(-1) ;
mpl_s = (1 - alpha_s)  y_sl_s ;
%
beta_s  mpk_s(1)  (ratio(1)  ps(1) + 1 - ratio(1)) + beta_s  lamd_s1(1)  (1 - delta) - lamd_s1 + lamd_s2  m_s q_s(1)  infl(1) R_s = 0 ; 
lamd_s1 = q_s  ratio  ;
% 
(ratio  ps + (1 - ratio) )  mpl_s = ratio  w_s ;
% 
ratio - beta_s  ratio(1)  R_s  infl(1) - lamd_s2 = 0 ;

%%%% 
c_b + b_s + b_e  + R_d(-1)  d(-1)  infl = b_s(-1)  R_s(-1)infl + b_e(-1)  R_e(-1)infl + d ;
d = m_b  (b_s + b_e) ;
lamd_b = c_b^(-sigma) ;
lamd_b - beta_b  lamd_b(1)  R_d infl(1) - llamd_b = 0 ;
-lamd_b + beta_b  lamd_b(1)  R_s infl(1) + m_bllamd_b = 0 ;
-lamd_b + beta_b  lamd_b(1)  R_e infl(1) + m_bllamd_b = 0 ;
 

%%%% 
i_e = ie - (fy_e  2)  (ieie(-1) - 1)^2  ie ; 
i_s = is - (fy_s  2)  (isis(-1) - 1)^2  is ;
1 = q_e  (1 - (fy_e  2)  (ieie(-1) - 1)^2   -  ie  fy_e  (ieie(-1) - 1) ie(-1) ) + sdf(1)  q_e(1) (fy_e  (ie(1)ie - 1) ) (ie(1)  ie)^2 ;
1 = q_s  (1 - (fy_s  2)  (isis(-1) - 1)^2   -  is  fy_s  (isis(-1) - 1) is(-1) ) + sdf(1)  q_s(1) (fy_s  (is(1)is - 1) ) (is(1)  is)^2 ;

%%%% 
z1 = y  x + thetesdf(+1)infl(+1)^(epsi)z1(+1);
z2 = y + thetesdf(+1)infl(+1)^(epsi-1)z2(+1);
inflp = (epsi(epsi-1))z1z2;
1 = theteinfl^(epsi-1) + (1-thete)inflp^(1-epsi);
sdf = beta_h  c(-1)c ;
vp = (1-thete)inflp^(-epsi) + theteinfl^(epsi)vp(-1);
y vp = y_e ; 

%%%% 政策
R_d Rs=(R_d(-1)Rs)^(rho_r)  (infl^(fy_pi)  (yy(-1))^(fy_y))^(1 - rho_r)  exp(mp) ;

%%%% 一般均衡
y = c + c_e + c_b + ie + is ;

A_sA_ss=(A_s(-1)A_ss)^(0.9)  exp(u_z + u_zs);
A_eA_es = (A_e(-1)A_es)^(0.9)  exp(u_z + u_ze);
mp = 0.9mp(-1) + u_mp ;

lny = log(y) ;
end ;

initval;
c =cs ;
l_s = l_ss ;
l_n = l_ns ;
d = ds;
c_e = c_es ;
i_e = i_es ;
b_e = b_es ;
k_e = k_es ;
y_e = y_es;
A_e = A_es ;
lamd_e = lamd_es ;
llamd_e = llamd_es ;
x = xs ;
mpl_e = mpl_es ;
mpk_e = mpk_es ; 
ratio = ratios ; 
xx = 0 ; 
i_s = i_ss;
b_s = b_ss ;
k_s = k_ss ;
y_s =y_ss ;
ps =pss ;
A_s = A_ss ;
lamd_s1 = lamd_s1s ;
lamd_s2 = lamd_s2s ;
mpl_s = mpl_ss ;
mpk_s = mpk_ss ;
c_b = c_bs ;
lamd_b = lamd_bs;
llamd_b = llamd_bs ;
is = iss ;
ie = ies ; 
z1 = z1s ;
z2 = z2s ; 
inflp = inflps ;
infl = infls;
vp = vps ;
y = ys ;
sdf =  sdfs ;
R_d = R_ds ;
R_e = R_es ;
R_s = R_ss ;
w_s = w_ss ;
w_n = w_ns ;
q_e = q_es;
q_s = q_ss; 
mp = 0 ; 
lny = 0 ;
end;

resid;
steady ;
check ;


shocks; 

var u_z = 0.0001 ;

end;



stoch_simul(periods = 1000,order=1,irf=50);