/*  garbaf(2013)                                                         */

var a bcb ce ch d etabar g ke kh lambdaco lambdae lambdah lcb lin 
    m me mh phie nwco nwin qk r ra rcb rd rin sigmaeta te xe xh ye yh y
    theta psicb la lce lch ld lke lkh lmh llcb llin lnwin lqk 
    lra lrcb lrd lrin lye lyh lambdaa lambdab lambdac lambdad;

varexo  eps_e eps_h eps_sigma;

parameters ALPHAE ALPHAH BETAE BETAH DELTA G jH jE K LTVE MU RHOE RHOH 
    RHOSIGMA STDERRA sigmaeta_ss V XI XICB w rd_ss a_ss bcb_ss ce_ss ch_ss d_ss
    etabar_ss g_ss ke_ss kh_ss lambdaco_ss lambdae_ss lambdah_ss 
    me_ss mh_ss m_ss lcb_ss lin_ss nwco_ss nwin_ss phie_ss psicb_ss 
    qk_ss r_ss ra_ss rcb_ss rin_ss te_ss theta_ss ye_ss yh_ss
    lambdaa_ss lambdab_ss lambdac_ss lambdad_ss PHI_ss frac; 

load chapt3ss.mat;
set_param_value('w',w)                  
set_param_value('ALPHAE',ALPHAE)        
set_param_value('ALPHAH',ALPHAH)        
set_param_value('BETAE',BETAE)          
set_param_value('BETAH',BETAH)          
set_param_value('DELTA',DELTA)          
set_param_value('jH',jH)                
set_param_value('jE',jE)                
set_param_value('K',K)                  
set_param_value('LTVE',LTVE)            
set_param_value('MU',MU)                
set_param_value('RHOE',RHOE)            
set_param_value('RHOH',RHOH)            
set_param_value('RHOSIGMA',RHOSIGMA)    
set_param_value('XI',XI)                
set_param_value('XICB',XICB)            
set_param_value('G',G)                  
set_param_value('STDERRA',STDERRA)      
set_param_value('V',V)                  
set_param_value('sigmaeta_ss',sigmaeta_ss)
set_param_value('rd_ss',rd_ss)
set_param_value('qk_ss',qk_ss) 
set_param_value('ke_ss',ke_ss) 
set_param_value('kh_ss',kh_ss) 
set_param_value('a_ss',a_ss) 
set_param_value('lin_ss',lin_ss) 
set_param_value('lcb_ss',lcb_ss)
set_param_value('rcb_ss',rcb_ss) 
set_param_value('etabar_ss',etabar_ss) 
set_param_value('ra_ss',ra_ss) 
set_param_value('sigmaeta_ss',sigmaeta_ss) 
set_param_value('ce_ss',ce_ss) 
set_param_value('ch_ss',ch_ss) 
set_param_value('mh_ss',mh_ss) 
set_param_value('me_ss',me_ss)
set_param_value('bcb_ss',bcb_ss) 
set_param_value('psicb_ss',psicb_ss) 
set_param_value('te_ss',te_ss) 
set_param_value('d_ss',d_ss) 
set_param_value('theta_ss',theta_ss) 
set_param_value('rin_ss',rin_ss) 
set_param_value('ye_ss',ye_ss)
set_param_value('yh_ss',yh_ss)
set_param_value('lambdae_ss',lambdae_ss) 
set_param_value('lambdah_ss',lambdah_ss) 
set_param_value('lambdaa_ss',lambdaa_ss)
set_param_value('lambdab_ss',lambdab_ss)
set_param_value('lambdac_ss',lambdac_ss)
set_param_value('lambdad_ss',lambdad_ss)
set_param_value('lambdaco_ss',lambdaco_ss) 
set_param_value('phie_ss',phie_ss) 
set_param_value('nwin_ss',nwin_ss)
set_param_value('g_ss',g_ss)
set_param_value('r_ss',r_ss)
%set_param_value('th_ss',th_ss)
set_param_value('nwco_ss',nwco_ss)
set_param_value('m_ss',m_ss)
set_param_value('PHI_ss',PHI_ss)
set_param_value('frac',frac) 

model;
%% Shocks
log(xe) = RHOE*log(xe(-1))+eps_e;   //Entrepreneurs Productivity 1
log(xh) = RHOH*log(xh(-1))+eps_h;   //Households Productivity 2
log(sigmaeta) = (1-RHOSIGMA)*log(sigmaeta_ss) + RHOSIGMA*log(sigmaeta(-1))+sigmaeta_ss*eps_sigma; //Standard deviation shock
%%
ch+d+qk*kh+mh+te=(1-DELTA)*qk*kh(-1)+mh(-1)+rd(-1)*d(-1)+yh;%+xh*(kh(-1))^(ALPHAH); 
(1/(ch)) = lambdah; 
lambdah = BETAH*lambdah(1)*rd; 
lambdah*qk = BETAH*(lambdah(1)*((1-DELTA)*qk(1)+ALPHAH*xh(1)*(kh)^(ALPHAH-1))); 
(jH/mh)=lambdah-BETAH*lambdah(1); 
%%
ce+ra(-1)*a(-1)+qk*ke+me+te=a+(1-DELTA)*qk*ke(-1)+me(-1)+ye;%+xe*(ke(-1))^(ALPHAE);
a = (((LTVE*(1-DELTA)*qk(1)*ke))/ra); 
lambdae-BETAE*lambdae(1)*ra = phie; 
(1/ce)=lambdae; 
lambdae*qk-BETAE*(lambdae(1)*((1-DELTA)*qk(1)+ALPHAE*xe(1)*(ke)^(ALPHAE-1)))=(((LTVE*(1-DELTA)*phie*qk(1)))/ra); 
(jE/me)=lambdae-BETAE*lambdae(1); 
%%
qk*K = qk*ke+qk*kh;
%%
nwco = rin(-1)*lin(-1) + theta(-1) - rd(-1)*d(-1);
lin + theta = d + nwco; 
theta = MU*d; 
lambdaco*(1-MU) = (rd-MU); 
lambdaco = rin; 
%%
nwin = frac*((1-(normcdf((((log(etabar)+((sigmaeta^2)/2))/(sigmaeta)))-sigmaeta)))*ra(-1)*a(-1) - (1-(normcdf((((log(etabar)+((sigmaeta^2)/2))/(sigmaeta))))))*(rcb(-1)*lcb(-1) + rin(-1)*lin(-1)));
a = lin + lcb + nwin;
etabar(1)*ra*a = rin*lin + rcb*lcb;
(1-(normcdf((((log(etabar(1))+((sigmaeta(1)^2)/2))/(sigmaeta(1)))))))*rcb*lcb + (1-XICB)*(1-w)*(normcdf((((log(etabar(1))+((sigmaeta(1)^2)/2))/(sigmaeta(1))))-sigmaeta(1)))*ra*a = r*lcb; 
(1-(normcdf((((log(etabar(1))+((sigmaeta(1)^2)/2))/(sigmaeta(1)))))))*rin*lin + (1-XI)*w*(normcdf((((log(etabar(1))+((sigmaeta(1)^2)/2))/(sigmaeta(1))))-sigmaeta(1)))*ra*a = rd*lin; 
%%
(1-(normcdf((((log(etabar(1))+((sigmaeta(1)^2)/2))/(sigmaeta(1))))-sigmaeta(1))))*ra - lambdaa + lambdab*ra*etabar(1) + lambdac*(1-XI)*(1-w)*(normcdf((((log(etabar(1))+((sigmaeta(1)^2)/2))/(sigmaeta(1))))-sigmaeta(1)))*ra + lambdad*(1-XI)*w*(normcdf((((log(etabar(1))+((sigmaeta(1)^2)/2))/(sigmaeta(1))))-sigmaeta(1)))*ra = 0;
lambdab*ra*a + lambdac*((etabar(1)*(normpdf((log(etabar(1))+(sigmaeta(1)^2)/2)/sigmaeta(1))/etabar(1)/sigmaeta(1)))*(1-XICB)*(1-w)*ra*a - (normpdf((log(etabar(1))+(sigmaeta(1)^2)/2)/sigmaeta(1))/etabar(1)/sigmaeta(1))*rcb*lcb) + lambdad*((etabar(1)*(normpdf((log(etabar(1))+(sigmaeta(1)^2)/2)/sigmaeta(1))/etabar(1)/sigmaeta(1)))*(1-XI)*w*ra*a - (normpdf((log(etabar(1))+(sigmaeta(1)^2)/2)/sigmaeta(1))/etabar(1)/sigmaeta(1))*rin*lin) = 0;
-(1-(normcdf((((log(etabar(1))+((sigmaeta(1)^2)/2))/(sigmaeta(1)))))))*rcb + lambdaa - lambdab*rcb + lambdac*((1-(normcdf((((log(etabar(1))+((sigmaeta(1)^2)/2))/(sigmaeta(1)))))))*rcb - r) = 0;
-(1-(normcdf((((log(etabar(1))+((sigmaeta(1)^2)/2))/(sigmaeta(1)))))))*rin + lambdaa - lambdab*rin + lambdad*((1-(normcdf((((log(etabar(1))+((sigmaeta(1)^2)/2))/(sigmaeta(1)))))))*rin - rd) = 0;
%%
m = me + mh; 
lcb + bcb = theta + m; 
psicb = r(-1)*bcb(-1) + rcb(-1)*lcb(-1) - theta(-1) - m(-1);
r = rd; 
rcb = rin;
%%
g  = 2*te + bcb - r(-1)*bcb(-1) + psicb;
g  = G*(ce+ch);
%%
ye = xe*(ke(-1))^(ALPHAE);
yh = xh*(kh(-1))^(ALPHAH);
y  = ye + yh;
%% Variables in LOGS
la = log(a)/STDERRA; 
lce = log(ce)/STDERRA; 
lch = log(ch)/STDERRA; 
ld  = log(d)/STDERRA;
lke  = log(ke)/STDERRA;
lkh  = log(kh)/STDERRA;
lmh = log(mh)/STDERRA;
llcb  = log(lcb)/STDERRA;
llin = log(lin)/STDERRA;
lnwin  = log(nwin)/STDERRA;
lqk  = log(qk)/STDERRA;
lra  = log(ra)/STDERRA;
lrcb  = log(rcb)/STDERRA;
lrd  = log(rd)/STDERRA;
lrin  = log(rin)/STDERRA;
lye  = log(ye)/STDERRA;
lyh  = log(yh)/STDERRA;

end;

initval;
xe = 1;
xh = 1;
sigmaeta = sigmaeta_ss;

a = a_ss; 
bcb = bcb_ss;
ce = ce_ss; 
ch = ch_ss; 
d  = d_ss;
etabar  = etabar_ss;
g  = g_ss;
ke  = ke_ss;
kh  = kh_ss;
lambdaa  = lambdaa_ss;
lambdab  = lambdab_ss;
lambdac  = lambdac_ss;
lambdad  = lambdad_ss;
lambdaco  = lambdaco_ss;
lambdae  = lambdae_ss;
lambdah  = lambdah_ss;
lcb  = lcb_ss;
lin = lin_ss;
me = me_ss;
mh = mh_ss;
m = m_ss;
phie  = phie_ss;
nwco = 0;
nwin  = nwin_ss;
qk  = qk_ss;
r  = r_ss;
ra  = ra_ss;
rcb  = rcb_ss;
rd  = rd_ss;
rin  = rin_ss;
te  = te_ss;
theta = theta_ss;
ye  = ye_ss;
yh  = yh_ss;
y = ye_ss+yh_ss;
psicb  = psicb_ss;

la = log(a_ss)/STDERRA; 
lce = log(ce_ss)/STDERRA; 
lch = log(ch_ss)/STDERRA; 
ld  = log(d_ss)/STDERRA;
lke  = log(ke_ss)/STDERRA;
lkh  = log(kh_ss)/STDERRA;
lmh = log(mh_ss)/STDERRA;
llcb  = log(lcb_ss)/STDERRA;
llin = log(lin_ss)/STDERRA;
lnwin  = log(nwin_ss)/STDERRA;
lqk  = log(qk_ss)/STDERRA;
lra = log(ra_ss)/STDERRA;
lrcb  = log(rcb_ss)/STDERRA;
lrd  = log(rd_ss)/STDERRA;
lrin  = log(rin_ss)/STDERRA;
lye  = log(ye_ss)/STDERRA;
lyh  = log(yh_ss)/STDERRA;
end;

resid(1);
steady;
check;

shocks;
    var eps_e ; stderr STDERRA;         % first
    var eps_h  ; stderr STDERRA;      % second
    var eps_sigma  ; stderr STDERRA;  % third
end;

stoch_simul(order=2,irf=50);