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 = 1/beta_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_e /k_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_e / infl(1) - llamd_e = 0 ;
 
%%%% 
ratio = exp(xx)/ (1 + exp(xx)) ;
xx= control_xx * log(y/y(-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_s/l_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)  *(ie/ie(-1) - 1)^2   -  ie * fy_e  (ie/ie(-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) * (is/is(-1) - 1)^2   -  is  *fy_s  (is/is(-1) - 1) /is(-1) ) + sdf(1) * q_s(1) *(fy_s  *(is(1)/is - 1) ) *(is(1)  /is)^2 ;

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

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

y = c + c_e + c_b + ie + is ;

A_s/A_ss=(A_s(-1)/A_ss)^(0.9) * exp(u_z + u_zs);
A_e/A_es = (A_e(-1)/A_es)^(0.9) * exp(u_z + u_ze);
mp = 0.9*mp(-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);