// Liu_wang_Zha
// This mod file is meant to add a constrained household to LWZ(2010): Do Credit Constrainys Amplify Macro Fluctuations ? 



//Define Variables 

var  muc mucb muh w ql R mue mub N Nc I Y Ch Ce Cc qk Lh Le Lc K B Bc ggamma gz gq lambdaz nuz lambdaq nuq lambdaa varphi psi theta   Nh S;


// Define Shocks 

varexo epsz epsnuz epsq epsnuq epsa epsvarphi epspsi epstheta ; 


// Define parameters 
parameters ggammap Omegahp Omegaep lambdaabarp lambdaqbarp gammahp gammaep betap Omegap lambdakp  thetabarp alphap phip deltap  rhozp rhoap rhonuzp rhoqp rhonuqp  rhovarphip rhopsip rhothetap R_ss mub_mue_ss ql_Le_Y_ss I_K_ss K_Y_ss Y_K_ss I_Y_ss B_Y_ss Ce_Y_ss Ch_Y_ss ql_Lh_Ch_ss Lh_Le_ss N_ss Le_Lbar_ss Lh_Lbar_ss Lc_Lbar_ss ql_Le_B_ss sigmaap sigmazp sigmanuzp sigmaqp sigmanuqp sigmavarphip sigmapsip sigmathetap betacp gammacp  phicp mucb_muc_ss Cc_Y_ss Bc_Y_ss w_Nc_Y_ss ql_Lc_Y_ss B_S_ss Bc_S_ss S_Y_ss ql_Lh_Y_ss w_Nh_Y_ss; 
           


// Load Parameters and steady states 
load LWZ_params.mat; 
load LWZ_ss.mat ; 



model;

0=-(ggammap-betacp*(1+lambdaabarp)*gammacp)*(ggammap-gammacp)*muc-(ggammap^2+gammacp^2*betacp*(1+lambdaabarp))*Cc+ggammap*gammacp*(Cc(-1)-ggamma)-betacp*lambdaabarp*gammacp*(ggammap-gammacp)*lambdaa(+1)+betacp*(1+lambdaabarp)*ggammap*gammacp*(Cc(+1)+ggamma(+1));
0=-w-muc+psi;
0=-(muc+ql)+betacp*phicp*mucb_muc_ss*ggammap*(mucb+ggamma(+1)+ql(+1))+betacp*(1+lambdaabarp)*(muc(+1)+ql(+1))+betacp*lambdaabarp*lambdaa(+1)+(1-phicp*betacp*(1+lambdaabarp)*(betap-betacp+1/phicp))*(varphi-Lc);
0=R/R_ss+(betacp/ggammap)*(1+lambdaabarp)*(muc(+1)-muc-ggamma(+1))+(betacp/ggammap)*lambdaabarp*lambdaa(+1)+mucb_muc_ss*(mucb-muc);

0= -Omegahp*muh-(ggammap^2+gammahp^2*betap*(1+lambdaabarp))*Ch + ggammap*gammahp*(Ch(-1)-ggamma)-betap*lambdaabarp*gammahp*(ggammap-gammahp)*lambdaa(+1)+betap*(1+lambdaabarp)*ggammap*gammahp*(Ch(+1)+ggamma(+1)); 
0=-w-muh+psi; 
0=-ql-muh+betap*(1+lambdaabarp)*(muh(+1)+ql(+1))+(1-betap*(1+lambdaabarp))*(varphi-Lh)+betap*lambdaabarp*lambdaa(+1); 
0=-muh+R+(muh(+1)+lambdaa(+1)*(lambdaabarp)/(1+lambdaabarp)-ggamma(+1));
0=-Omegaep*mue-(ggammap^2+betap*gammaep^2)*Ce+ggammap*gammaep*(Ce(-1)-ggamma)+betap*ggammap*gammaep*(Ce(+1)+ggamma(+1)); 

0=-w+Y-N;
0=-qk+(1+betap)*Omegap*lambdakp^2*I-Omegap*lambdakp^2*I(-1)+Omegap*lambdakp^2*(ggamma+gq)-betap*Omegap*lambdakp^2*(I(+1)+ggamma(+1)+gq(+1)); 

0=-qk-mue+((mub_mue_ss*thetabarp)/(lambdaqbarp))*(mub+theta)+(betap*(1-deltap)/lambdakp)*(qk(+1)-gq(+1)-ggamma(+1))+(1-(mub_mue_ss*thetabarp/(lambdaqbarp)))*mue(+1)+((mub_mue_ss*thetabarp)/(lambdaqbarp))*(qk(+1)-gq(+1))+betap*alphap*(1-phip)*(Y_K_ss)*(Y(+1)-K); 

0=-ql-mue+(mub_mue_ss)*ggammap*thetabarp*(theta+mub)+(1-mub_mue_ss*ggammap*thetabarp)*mue(+1)+mub_mue_ss*ggammap*thetabarp*(ql(+1)+ggamma(+1))+betap*ql(+1)+(1-betap-betap*lambdaabarp*thetabarp)*(Y(+1)-Le) ;

0=-mue+R+(1/(1+lambdaabarp))*((mue(+1)-ggamma(+1))+lambdaabarp*mub);

0=-Y+alphap*phip*Le(-1)+alphap*(1-phip)*K(-1)+(1-alphap)*N- ((1-phip)*alphap/(1-(1-phip)*alphap))*(gz+gq);

0=-K+ (1-deltap)/(lambdakp)*(K(-1)-ggamma-gq)+(1-((1-deltap)/(lambdakp)))*I; 

0=-B+theta+ggammap*thetabarp*(ql_Le_B_ss)*(ql(+1)+Le+ggamma(+1))+(1-ggammap*thetabarp*ql_Le_B_ss)*(qk(+1)+K-gq(+1));

0=-alphap*Y+Ce*Ce_Y_ss+I*I_Y_ss+(Le-Le(-1))*ql_Le_Y_ss+((1*B_Y_ss)/(ggammap))*(B(-1)-ggamma)-((1*B_Y_ss)/(R_ss))*(B-R);

0=Cc_Y_ss*Cc+ql_Lc_Y_ss*(Lc-Lc(-1))-(Bc_Y_ss/R_ss)*(Bc-R)-w_Nc_Y_ss*(w+Nc)+(Bc_Y_ss/ggammap)*(Bc(-1)-ggamma);

0= Ch_Y_ss*Ch+ql_Lh_Y_ss*(Lh-Lh(-1))-w_Nh_Y_ss*(w+Nh)-(S_Y_ss/ggammap)*(S(-1)-ggamma)+(S_Y_ss/R_ss)*(S-R);

0=-(1-alphap)*N+w_Nc_Y_ss*Nc+w_Nh_Y_ss*Nh; 

0=-Bc+ql(+1)+Lc+ggamma(+1); 

0=Lh*Lh_Lbar_ss+Lc*Lc_Lbar_ss+Le*Le_Lbar_ss;
0=-S+B_S_ss*B+Bc_S_ss*Bc; 

// 0=-Y+Ch*Ch_Y_ss+ Cc*Cc_Y_ss+Ce*Ce_Y_ss+I*I_Y_ss; 






0=-gz+lambdaz+nuz-nuz(-1); 
0=-gq+lambdaq+nuq-nuq(-1); 
0=-ggamma+1/(1-(1-phip)*alphap)*gz+(1-phip)*alphap/(1-(1-phip)*alphap)*gq;
 
0=-lambdaz+rhozp*lambdaz(-1)+epsz;
0=-nuz+rhonuzp*nuz(-1)+epsnuz;
0=-lambdaq+rhoqp*lambdaq(-1)+epsq;
0=-nuq+rhonuqp*nuq(-1)+epsnuq;
0=-lambdaa+rhoap*lambdaa(-1)+epsa;
0=-varphi+rhovarphip*varphi(-1)+epsvarphi;
0=-psi+rhopsip*psi(-1)+epspsi;
0=-theta+rhothetap*theta(-1)+epstheta; 
end; 

shocks;
 var epsz; stderr sigmazp;
//  var epsnuz; stderr sigmanuzp;
 // var epsq;  stderr sigmaqp;
 var epsnuq; stderr sigmanuqp;
// var epsa; stderr sigmaap;
// var epsc; stderr sigmaap; 

var epsvarphi; stderr sigmavarphip ;
// var epspsi; stderr sigmapsip;
 var epstheta; stderr sigmathetap;
end;

initval;
muh =0;
w  =0;
ql  =0;
R =0;
mue =0;
mub  =0;
N =0;
I =0;
Y =0;
Ch =0;
Ce =0;
qk =0;
Lh  =0;
Le =0;
K  =0;
B =0;
ggamma =0;
gz =0;
gq  =0;
lambdaz =0;
nuz =0;
lambdaq =0;
nuq =0;
lambdaa =0;
varphi =0; 
psi  =0;
theta =0;
Cc=0;
Nh=0;
Nc=0;
S=0;
Bc=0;

end;


steady;


check;

stoch_simul(order = 1);