% Steady state for OPEN ECONOMY 
% No Mpru policy

clear all
clc

format compact
format long

% First, we declare the parameters we have directly from the papers
betta = 0.99   %Discount factor
hcons = 0.5 %Habit persistence 0 in GS and 0.9 in JFV
vartheta = 0.5  %Frisch elasticity of labor4
vareps= 8.577 %10      %Elasticity of substitution across intermediate goods
theta = 0.8     %Calvo pricing parameter
thetafor = 0.8  %Calvo pricing parameter FOREIGN
chi   = 0.6  %Degree of indexation
alfa  = 0.33%0.22 %   %Capital share intermediate production
delta = 0.023%0.01%    %Capital depreciation rate

pdef  = 0.03    %Annual probability of default (BGG)
Fwbar = 1-((1-pdef)^0.25) %Implied quarterly probability of default
pdeffor  = 0.03    %Annual probability of default (BGG) FOREIGN
Fwbarfor = 1-((1-pdeffor)^0.25) %Implied quarterly probability of default FOREIGN

Pihss=1
Pifss=1

muc   = 0.15%2    %bankruptcy costs
mucfor   = 0.15%2    %bankruptcy costs FOREIGN

%Fiscal variables
g_y   = 0.2%Spain in 2007 was 0.39. This ensures a reasonable tauc
    taul  = 0.24%0.24    %Labor income tax for Spain 2007
    taur  = 0.42%0.32    %Capital income tax for Spain 2007
gfor_yfor = 0.2%0.3 %FOREIGN France in 2007 was 0.52 in g_y but with higher d_y. This ensures a reasonable taucfor
    taulfor  = 0.24%0.32%Labor income tax for France 2007
    taurfor  = 0.42%0.42  %Capital income tax for France 2007

%Measured steady states from data
Piss  = 1   %Target gross inflation
 
labss = 1/3     %Time devoted to work
labssfor = 1/3     %Time devoted to work FOREIGN
qss   = 1       %Tobin's q
qssfor   = 1       %Tobin's q FOREIGN

%Eqn #3
Rdss   = Piss/betta  %Gross nominal interest rate
%Rdss  = Rss     %Gross nominal interest rate on deposits when taur = 0
  
%Eqn #2
Rss    = ((Rdss - 1)/(1-taur)) + 1
%Eqn #35
Rssfor = Rss %FOREIGN
%Eqn #57
Rdssfor = (Rssfor-1)*(1-taurfor)+1 %Gross nominal interest rate FOREIGN
%Eqn #58
Pissfor = betta*Rdssfor
omegaB =  0.0043 %Quint and Rabanal %0.000742%Schmitt-Grohé and Uribe (2003)%0%complete markets%
checkvec(1) = Rss - Rdss;
checkvec(2) = Rss - Rssfor;
checkvec(3) = Rssfor - Rss;

%Goods market open economy parameters
size = 0.5 %Size of the home country with respect to the rest of th union
zet = 1.5 % Elasticity of substitution between domestic and foreign goods. Ester Faia 1.5.
%(1-varphi) is the home bias of the home country and varphifor is the home
%bias of the foreign country

%Using n/(k-n)=2 and eqn#14 we get
bbar_k = 1/3  %0.35 %Recall bbar = b/p
bbar_kfor= 1/3 %1/4 %FOREIGN

n_k    = 1- bbar_k 
n_kfor    = 1- bbar_kfor %FOREIGN 

%Eqn #18 and Eqn#70
Pistar = ((1-theta*(Pihss^((vareps-1)*(1-chi))))/(1-theta))^(1/(1-vareps))
Pistarfor = ((1-thetafor*(Pifss^((vareps-1)*(1-chi))))/(1-thetafor))^(1/(1-vareps)) %FOREIGN

%Eqns #15-17 and #67-69
mcss = ((vareps-1)/vareps)*(Pistar*(1-betta*theta*(Pihss^(vareps*(1-chi)))))/((1-betta*theta*(Pihss^((vareps-1)*(1-chi)))))
mcssfor = ((vareps-1)/vareps)*(Pistarfor*(1-betta*thetafor*(Pifss^(vareps*(1-chi)))))/((1-betta*thetafor*(Pifss^((vareps-1)*(1-chi)))))

%Eqn #41 and #90
nuss   = ((1-theta)*(Pistar^(-vareps)))/(1-theta*(Pihss^(vareps*(1-chi)))) 
nussfor   = ((1-thetafor)*(Pistarfor^(-vareps)))/(1-thetafor*(Pifss^(vareps*(1-chi)))) 

pars = [Fwbar;muc;bbar_k];
x0 = [0.5;0.5];
options = optimset('TolX',1e-14);
[x0res,fvalw] = fminsearch('fomega2',x0,options,pars);
omegabar = x0res(1)
sigmaw = x0res(2)
fvalw;
muw   = -(1/2)*(sigmaw^2)
Fwbar = logncdf(omegabar,muw,sigmaw)
Gwbar = 1 - normcdf(((-log(omegabar) + 0.5*(sigmaw^2))/sigmaw))
Fwbar_w = lognpdf(omegabar,muw,sigmaw)

%FOREIGN
parsfor = [Fwbarfor;mucfor;bbar_kfor];
x0for = [0.5;0.5];
options = optimset('TolX',1e-14);
[x0forres,fvalwfor] = fminsearch('fomega2open',x0for,options,parsfor);
omegabarfor = x0forres(1)
sigmawfor = x0forres(2)
fvalwfor;
muwfor   = -(1/2)*(sigmawfor^2)
Fwbarfor = logncdf(omegabarfor,muwfor,sigmawfor)
Gwbarfor = 1 - normcdf(((-log(omegabarfor) + 0.5*(sigmawfor^2))/sigmawfor))
Fwbar_wfor = lognpdf(omegabarfor,muwfor,sigmawfor)
%FOREIGN

Gamawbar = omegabar*(1-Fwbar) + Gwbar
Gamawbarfor = omegabarfor*(1-Fwbarfor) + Gwbarfor %FOREIGN

Gamawbar_w = 1-Fwbar
Gamawbar_wfor = 1-Fwbarfor %FOREIGN
Fwbar_w2 = (1/(omegabar*sigmaw*((2*pi)^0.5)))*exp(-((log(omegabar) - muw)^2)/(2*(sigmaw^2)))
Fwbar_w2for = (1/(omegabarfor*sigmawfor*((2*pi)^0.5)))*exp(-((log(omegabarfor) - muwfor)^2)/(2*(sigmawfor^2))) %FOREIGN
Gwbar_w = omegabar*Fwbar_w
Gwbar_wfor = omegabarfor*Fwbar_wfor %FOREIGN

Fwbar_ww = -(Fwbar_w/omegabar)*(1+((log(omegabar) - muw)/(sigmaw^2)))
Fwbar_wsig = (1/(omegabar*(sigmaw^2)*((2*pi)^0.5)))*(exp(-((log(omegabar) - muw)^2)/(2*(sigmaw^2))))*((((4*(log(omegabar))^2) - sigmaw^4)/((4*(sigmaw^2)))) - 1)
Fwbar_wwfor = -(Fwbar_wfor/omegabarfor)*(1+((log(omegabarfor) - muwfor)/(sigmawfor^2))) %FOREIGN
Fwbar_wsigfor = (1/(omegabarfor*(sigmawfor^2)*((2*pi)^0.5)))*(exp(-((log(omegabarfor) - muwfor)^2)/(2*(sigmawfor^2))))*((((4*(log(omegabarfor))^2) - sigmawfor^4)/((4*(sigmawfor^2)))) - 1) %FOREIGN

Fwbar_sig = normpdf((log(omegabar) + 0.5*(sigmaw^2))/sigmaw)*(1 - ((log(omegabar) + 0.5*(sigmaw^2))/(sigmaw^2)));
Gwbar_sig = normpdf(((log(omegabar) + 0.5*(sigmaw^2))/sigmaw) - sigmaw)*(1 - ((log(omegabar) + 0.5*(sigmaw^2))/(sigmaw^2)) - 1);
Fwbar_sigfor = normpdf((log(omegabarfor) + 0.5*(sigmawfor^2))/sigmawfor)*(1 - ((log(omegabarfor) + 0.5*(sigmawfor^2))/(sigmawfor^2))); %FOREIGN
Gwbar_sigfor = normpdf(((log(omegabarfor) + 0.5*(sigmawfor^2))/sigmawfor) - sigmawfor)*(1 - ((log(omegabarfor) + 0.5*(sigmawfor^2))/(sigmawfor^2)) - 1); %FOREIGN

Gamawbar_sig = -omegabar*Fwbar_sig + Gwbar_sig;
Gamawbar_sigfor = -omegabarfor*Fwbar_sigfor + Gwbar_sigfor; %FOREIGN

% Additional variables required in the log-linearization
Psi1 = 1-Gamawbar %ok
Psi1w = -(1 - Fwbar) %ok
Psi1sigmaw = -Gamawbar_sig
Psi1wb = -Gamawbar_w; %ok
%disp('check')
checkvec(4) = Psi1w - Psi1wb;

Psi1for = 1-Gamawbarfor %FOREIGN
Psi1wfor = -(1 - Fwbarfor) %FOREIGN
Psi1sigmawfor = -Gamawbar_sigfor %FOREIGN
Psi1wbfor = -Gamawbar_wfor; %FOREIGN
%disp('check')
checkvec(5) = Psi1wfor - Psi1wbfor; %FOREIGN

Psi2 = (1-Fwbar)/(1-Fwbar - muc*omegabar*Fwbar_w) %p.27 JFV, en la p.28 es erroneo
Psi2w = (-Fwbar_w*(1-Fwbar-muc*omegabar*Fwbar_w)-(1-Fwbar)*(-Fwbar_ww*muc*omegabar - muc*Fwbar_w- Fwbar_w))/(1-Fwbar-muc*omegabar*Fwbar_w)^2
Psi2sigmaw = ((1-Fwbar)*muc*omegabar*Fwbar_wsig + Fwbar_sig*(1-Fwbar)-((1-Fwbar - muc*omegabar*Fwbar_w)*Fwbar_sig))/(1-Fwbar-muc*omegabar*Fwbar_w)^2
Psi2sigmaw2 = (-Fwbar_sig*(1-Fwbar-muc*omegabar*Fwbar_w)-(1-Fwbar)*(-Fwbar_sig-muc*omegabar*Fwbar_wsig))/(1-Fwbar-muc*omegabar*Fwbar_w)^2
Psi2wb = ((1-Fwbar)*Fwbar_w + muc*Fwbar_w + muc*omegabar*Fwbar_ww - Fwbar*muc*Fwbar_w- Fwbar*muc*omegabar*Fwbar_ww-(1-Fwbar - muc*omegabar*Fwbar_w)*Fwbar_w)/(1-Fwbar-muc*omegabar*Fwbar_w)^2;
checkvec(6)= Psi2w - Psi2wb;
checkvec(7)= Psi2sigmaw - Psi2sigmaw2;


%FOREIGN
Psi2for = (1-Fwbarfor)/(1-Fwbarfor - mucfor*omegabarfor*Fwbar_wfor) %p.27 JFV
Psi2wfor = (-Fwbar_wfor*(1-Fwbarfor-muc*omegabarfor*Fwbar_wfor)-(1-Fwbarfor)*(-Fwbar_wwfor*muc*omegabarfor - muc*Fwbar_wfor- Fwbar_wfor))/(1-Fwbarfor-muc*omegabarfor*Fwbar_wfor)^2
Psi2sigmawfor = ((1-Fwbarfor)*muc*omegabarfor*Fwbar_wsigfor + Fwbar_sigfor*(1-Fwbarfor)-((1-Fwbarfor - muc*omegabarfor*Fwbar_wfor)*Fwbar_sigfor))/(1-Fwbarfor-muc*omegabarfor*Fwbar_wfor)^2
Psi2sigmaw2for = (-Fwbar_sigfor*(1-Fwbarfor-muc*omegabarfor*Fwbar_wfor)-(1-Fwbarfor)*(-Fwbar_sigfor-muc*omegabarfor*Fwbar_wsigfor))/(1-Fwbarfor-muc*omegabarfor*Fwbar_wfor)^2
Psi2wbfor = ((1-Fwbarfor)*Fwbar_wfor + muc*Fwbar_wfor + muc*omegabarfor*Fwbar_wwfor - Fwbarfor*muc*Fwbar_wfor- Fwbarfor*muc*omegabarfor*Fwbar_wwfor-(1-Fwbarfor - muc*omegabarfor*Fwbar_wfor)*Fwbar_wfor)/(1-Fwbarfor-muc*omegabarfor*Fwbar_wfor)^2;
checkvec(8)= Psi2wfor - Psi2wbfor;
checkvec(9)= Psi2sigmawfor - Psi2sigmaw2for;
%FOREIGN

Psi3 = Gamawbar - muc*Gwbar
Psi3w = 1-Fwbar - muc*omegabar*Fwbar_w
Psi3sigmaw = Gamawbar_sig - muc*Gwbar_sig
Psi3for = Gamawbarfor - mucfor*Gwbarfor %FOREIGN
Psi3wfor = 1-Fwbarfor - mucfor*omegabarfor*Fwbar_wfor %FOREIGN
Psi3sigmawfor = Gamawbar_sigfor - mucfor*Gwbar_sigfor %FOREIGN

omega_a = ((Psi1w*omegabar)/Psi1) - ((Psi2w*omegabar)/Psi2)
sigma_a = ((Psi1sigmaw*sigmaw)/Psi1) - ((Psi2sigmaw*sigmaw)/Psi2)
omega_afor = ((Psi1wfor*omegabarfor)/Psi1for) - ((Psi2wfor*omegabarfor)/Psi2for) %FOREIGN
sigma_afor = ((Psi1sigmawfor*sigmawfor)/Psi1for) - ((Psi2sigmawfor*sigmawfor)/Psi2for); %FOREIGN

omega_b = ((Psi3w*omegabar)/Psi3)
sigma_b = ((Psi3sigmaw*sigmaw)/Psi3)
omega_bfor = ((Psi3wfor*omegabarfor)/Psi3for) %FOREIGN
sigma_bfor = ((Psi3sigmawfor*sigmawfor)/Psi3for) %FOREIGN

omega_c = (((omegabar^2)*Fwbar_w)/Gwbar)
sigma_c = ((Gwbar_sig*sigmaw)/Gwbar)
omega_cfor = (((omegabarfor^2)*Fwbar_wfor)/Gwbarfor) %FOREIGN
sigma_cfor = ((Gwbar_sigfor*sigmawfor)/Gwbarfor) %FOREIGN

%disp('check')
checkvec(10) = Fwbar - logncdf(omegabar,muw,sigmaw);

itass=1
itassfor=1

%Eqn. #26 and #78
Rkss = (bbar_k*Rss*itass)/(Gamawbar - muc*Gwbar)
Rkssfor = (bbar_kfor*Rssfor*itassfor)/(Gamawbarfor - mucfor*Gwbarfor)%FOREIGN

%With wbar back to get Rk from eqn 13 for example
checkvec(11) = Rkss - (Rss*bbar_k*itass)/(Gamawbar - muc*Gwbar);
%break
%Eqn #24 and #76
Rlss = omegabar*Rkss/bbar_k
Rlssfor = omegabarfor*Rkssfor/bbar_kfor

%Eqn #23 and #75
rss = (Rkss/Piss) -1 + delta
rssfor = (Rkssfor/Pissfor) -1 + delta;

%Eqn #14 and #66
wss = (mcss/(((1/(1-alfa))^(1-alfa))*((1/alfa)^(alfa))*(rss^alfa)))^(1/(1-alfa))
wssfor = (mcssfor/(((1/(1-alfa))^(1-alfa))*((1/alfa)^(alfa))*(rssfor^alfa)))^(1/(1-alfa)) %FOREIGN

%Eqn #13 and #65
kss = (alfa/(1-alfa))*(wss/rss)*labss
kssfor = (alfa/(1-alfa))*(wssfor/rssfor)*labssfor %FOREIGN

bbar = bbar_k*kss
bbarfor = bbar_kfor*kssfor %FOREIGN
nss = n_k*kss
nssfor = n_kfor*kssfor %FOREIGN

%Eqn #13 and #65
k_l = (alfa/(1-alfa))*(wss/rss)
k_lfor = (alfa/(1-alfa))*(wssfor/rssfor)
%disp('check')
checkvec(12) = wss - rss*((1-alfa)/alfa)*k_l;
%Eqn #19 and #71
iss = delta*kss
issfor = delta*kssfor
%Eqn #40 and #89
yss = (1/nuss)*(kss^alfa)*(labss^(1-alfa))
yssfor = (1/nussfor)*(kssfor^alfa)*(labssfor^(1-alfa))
%Eqn #40 and #89
k_y = nuss*(k_l^(1-alfa))
k_yfor = nussfor*(k_lfor^(1-alfa))
%Eqn #19 and #71
i_y = delta*k_y
ifor_yfor = delta*k_yfor
checkvec(13) = i_y - (iss/yss);
checkvec(14) = k_y - (kss/yss);
checkvec(15) = k_yfor - (kssfor/yssfor);

dss_y = 0.6
dss_yfor = 0.6 %FOREIGN
dss4 = dss_y*yss
dss4for = dss_yfor*yssfor %FOREIGN
dssEMU=(dss4^size)*dss4for^(1-size)
yssEMU=(yss^size)*yssfor^(1-size)
gss = g_y*yss
gssfor = gfor_yfor*yssfor %FOREIGN

tradess=1
%indexss^(1-zet) = (1-varphi)+varphi*(tradess)^(1-zet)
indexss=1
%indexforss^(1-zet) = (varphifor)+(1-varphifor)*(1/tradess)^(1-zet)
indexforss=1

chfor_yfor=0.32%0.11
cf_y=0.27%0.1%
cfss=cf_y*yss
chssfor=chfor_yfor*yssfor

netxss = (1-size)*chssfor-size*cfss
netxforss = size*cfss - (1-size)*chssfor

nx_y=(1-size)*chfor_yfor*yssfor/yss-size*tradess*cf_y
nxfor_yfor=size*cf_y*yss/yssfor - (1-size)*(1/tradess)*chfor_yfor

css=indexss*(yss-iss-gss-muc*kss*Gwbar*(rss+1-delta)-netxss/size)
cssfor=indexforss*(yssfor-issfor-gssfor-mucfor*kssfor*Gwbarfor*(rssfor+1-delta)-netxforss/(1-size))

varphi=cfss/(((indexss*tradess)^(-zet))*css)
varphifor=1-(chssfor/(((indexforss/tradess)^(-zet))*cssfor))

chss = (1-varphi)*css*(indexss)^(-zet)
cfssfor = varphifor*cssfor*(indexforss)^(-zet)

bondss= netxss/(size*(((Rss)/Piss)-1))%NFAS
bondforss=((-size)/(1-size))*bondss/tradess%zero net supply in real terms

c_y=css/yss
cfor_yfor=cssfor/yssfor
tauc = ((((Rdss/Piss)-1)*dss_y)+g_y - (taul*(wss*labss/yss) + taur*((Rss-1)/Pihss)*((bbar-bondss)/yss)))/c_y
taucfor = ((((Rdssfor/Pissfor)-1)*dss_yfor)+gfor_yfor - (taulfor*(wssfor*labssfor/yssfor) + taurfor*((Rssfor-1)/Pifss)*((bbarfor-bondforss)/yssfor)))/cfor_yfor %FOREIGN
    tax_y4 = (tauc*c_y + taul*(wss*labss/yss) + taur*((Rss-1)/Piss)*((bbar-bondss)/yss))
    tax_y4for = (taucfor*cfor_yfor + taulfor*(wssfor*labssfor/yssfor) + taurfor*((Rssfor-1)/Pissfor)*((bbarfor-bondforss)/yssfor))%FOREIGN
    taxss4 = tax_y4*yss
    taxss4for = tax_y4for*yssfor %FOREIGN
    gss4 = dss4*(1-(Rdss/Piss)) + taxss4
    gss4for = dss4for*(1-(Rdssfor/Pissfor)) + taxss4for %FOREIGN
    taxss = taxss4
    taxssfor = taxss4for %FOREIGN
    gcoef = gss4
    gcoeffor = gss4for %FOREIGN
    checkvec(16) = dss4 - (Rdss/Piss)*dss4 - gss4+taxss4

%back to eqn #1 to get lpsi
lpsi = ((1-betta*hcons)/((1-hcons)*css))*((1-taul)/(1+tauc))*(wss/(labss^vartheta))
lpsifor = ((1-betta*hcons)/((1-hcons)*cssfor))*((1-taulfor)/(1+taucfor))*(wssfor/(labssfor^vartheta))
w_r = wss/rss
w_rfor = wssfor/rssfor

ch_y=chss/yss
chfor_y=chssfor/yss

checkvec(17) = chss - (css/indexss-tradess*indexss*cfss)
checkvec(18) = cfssfor - (cssfor/indexforss-chssfor*indexforss/tradess)

%Eqn9
wss2 = ((mcss*(w_r^alfa))/(((1/(1-alfa))^(1-alfa))*((1/alfa)^alfa)))
wss2for = ((mcssfor*(w_rfor^alfa))/(((1/(1-alfa))^(1-alfa))*((1/alfa)^alfa)))
%disp('check')
checkvec(19) = wss - wss2;

%Eqn 23
%disp('check')
checkvec(20) = k_l - (k_y/nuss)^(1/(1-alfa));
%using l=1/3
kss2 = k_l*labss
kss2for = k_lfor*labssfor
checkvec(21) =kss - kss2;
checkvec(22) =kssfor - kss2for;

%disp('check')
%nuss - k_y/(((alfa/(1-alfa))*(wss/rss))^(1-alfa))
checkvec(23) = nuss - k_y/(k_l^(1-alfa));
checkvec(24) = k_l - w_r*(alfa/(1-alfa));
checkvec(25) = k_l - (wss/rss)*(alfa/(1-alfa));
%disp('check')
checkvec(26) = k_l - ((alfa/(1-alfa))*(wss/rss));
%break
%alternative: IT DOES NOT WORK, C_Y =2!!!!
%lpsi = 3.83
%c_y = ((1-betta*hcons)/(1-hcons))*((1-taul)/(1+tauc))*((1-alfa)/alfa)*(k_y*rss/(lpsi*(labss^(vartheta+1))))
%break

%Eqn #1 and #56
lambdass = ((1-betta*hcons)/(1-hcons))*(1/css)*(1/(1+tauc))
lambdassfor = ((1-betta*hcons)/(1-hcons))*(1/cssfor)*(1/(1+taucfor))

%Eqn #15 and #67
f1ss = (lambdass*mcss*yss)/(1-betta*theta*(Pihss^(vareps*(1-chi))))
f1ssfor = (lambdassfor*mcssfor*yssfor)/(1-betta*thetafor*(Pifss^(vareps*(1-chi))))
%Eqn #16 and #68
f2ss = (lambdass*Pistar*yss)/(1-betta*theta*(Pihss^((vareps-1)*(1-chi))))
f2ssfor = (lambdassfor*Pistarfor*yssfor)/(1-betta*thetafor*(Pifss^((vareps-1)*(1-chi))))
%Eqn8
wss3 = k_l*rss*((1-alfa)/alfa)
wss3for = k_lfor*rssfor*((1-alfa)/alfa)
checkvec(27) = wss - wss3;
checkvec(28) = wssfor - wss3for;

%Eqn #5
checkvec(29) = vareps*f1ss - (vareps-1)*f2ss;
checkvec(30) = vareps*f1ssfor - (vareps-1)*f2ssfor;

checkvec(33) = Fwbar_w - Fwbar_w2
checkvec(34) = Fwbar_wfor - Fwbar_w2for

checkvec(35) = yss -(chss+((chssfor)*(1-size)/size)+iss+gss+muc*Gwbar*(1+rss-delta)*kss)
checkvec(36) = 1 -(ch_y+(chfor_y)*(1-size)/size+i_y+g_y+muc*Gwbar*(1+rss-delta)*k_y)
checkvec(37) = yssfor-((size/(1-size))*(cfss)+cfssfor+issfor+gssfor+mucfor*kssfor*Gwbarfor*(rssfor+1-delta))
checkvec(38) = yssfor/yssfor-((size/(1-size))*(cfss/yssfor)+(cfssfor/yssfor)+gssfor/yssfor + (issfor/yssfor)+mucfor*(kssfor/yssfor)*Gwbarfor*(rssfor+1-delta))
checkvec(39) = css -( (((1-varphi)^(1/zet))*(chss^((zet-1)/zet))+(varphi^(1/zet))*(cfss^((zet-1)/zet)))^(zet/(zet-1)))
checkvec(40) = cssfor - ((((1-varphifor)^(1/zet))*(chssfor^((zet-1)/zet))+(varphifor^(1/zet))*(cfssfor^((zet-1)/zet)))^(zet/(zet-1)))
checkvec(41) = cfss-varphi*css*(tradess*indexss)^(-zet)
checkvec(42) = chssfor-(1-varphifor)*cssfor*(indexforss/tradess)^(-zet)

gamma_ebar = 3.67
gamma_e = 1/(1+exp(-gamma_ebar))
gamma_ebarfor = 3.67
gamma_efor = 1/(1+exp(-gamma_ebarfor))

a1 = (gamma_e/(Piss*nss))*((1-muc*Gwbar)*Rkss*kss - Rss*bbar*itass)
a2 = -(gamma_e/(Piss*nss))*muc*Gwbar*Rkss*kss
a3 = (gamma_e/(Piss*nss))*(1-muc*Gwbar)*Rkss*kss
a4 = -(gamma_e/(Piss*nss))*Rss*bbar*itass
a1for = (gamma_efor/(Pissfor*nssfor))*((1-mucfor*Gwbarfor)*Rkssfor*kssfor - Rssfor*bbarfor*itassfor)
a2for = -(gamma_efor/(Pissfor*nssfor))*mucfor*Gwbarfor*Rkssfor*kssfor
a3for = (gamma_efor/(Pissfor*nssfor))*(1-mucfor*Gwbarfor)*Rkssfor*kssfor
a4for = -(gamma_efor/(Pissfor*nssfor))*Rssfor*bbarfor*itassfor

wsse = nss - gamma_e*(1/Piss)*(Rkss*kss*(1-muc*Gwbar) - Rss*bbar)
if wsse<0;
    disp('WARNING: Negative wsse home');
    pause(2);
end;
wssefor = nssfor - gamma_efor*(1/Pissfor)*(Rkssfor*kssfor*(1-mucfor*Gwbarfor) - Rssfor*bbarfor)
if wssefor<0;
    disp('WARNING: Negative wsse foreign');
    pause(2);
end;

S2p1 = 14.477%4.75%
S2p1for = 14.477%4.75%

%Now it's time to choose the macroprudential scenario
disp('1. MpruRatio HOME')
disp('2. MpruRatio FOREIGN')
disp('3. MpruRatio HOME & FOREIGN')
disp('4. NoMpru')
mpru = input('Which mpru scenario do you want?')
if mpru==1;
    gamma_ita=1.75
    gamma_itafor=0
elseif mpru==2;
    gamma_ita=0
    gamma_itafor=1.75
elseif mpru==3;
    gamma_ita=1.75
    gamma_itafor=1.75
elseif mpru==4;
    gamma_ita=0
    gamma_itafor=0
end;


disp('checks')
checkvec'

%Parameters for policy rules/shock processes
%Gov't expenditure
gamma_g = 0.95;
%d_g  = -0.09;
%d_y  = 0;
%sigma_g  =;

    %tax on consumption
    gamma_c  = 0.95;
    %sigma_tauc  =;
    %Tax on labor income
    gamma_l  = 0.95;
    %sigma_taul =;
    %Tax on deposit returns
    gamma_r  = 0.95;
    %sigma_taur  =;

%Intertemporal shock
rhod  = 0.95;
%sigma_phi  =;
%Technology shocks
rhoa  = 0.95;
%sigma_z  =;
%Dispersion process
rhosigma  = 0.95;
%eta_sigma =;

%sigma_s  =;
%Entrepreneurs entry
rhogammae  = 0.95;
%sigma_gammae =;

%Taylor rule
%gamma_R = 0.95;
%gamma_pi = 1.5;%0.95;
%gamma_y = 0%.95;4

%Gov't budget constraint1

taxcoef = taxss 
taxcoeffor = taxssfor %FOREIGN

%Now save the main values for .mod
if mpru==1
save NI_paramfileopenMPRUHOME betta hcons vartheta theta vareps chi alfa delta Piss Rss Pistar...
    S2p1 Rkss omega_a sigma_a omega_b sigma_b qss rss nss kss bbar_k... 
    a1 gamma_e a2 omega_c sigma_c a3 a4 gcoef taxcoef tauc taul taur css wss labss bbar c_y i_y g_y dss_y yss...
    muc Gwbar gamma_g gamma_c gamma_l gamma_r rhod rhoa rhosigma rhogammae dss4 Rdss Rlss...
    omegabar Fwbar_w Fwbar_sig Gwbar_w Gwbar_sig sigmaw zet varphi varphifor omegaB...
    thetafor Pissfor Rssfor Pistarfor S2p1for Rkssfor omega_afor sigma_afor...
    omega_bfor sigma_bfor qssfor rssfor nssfor kssfor bbar_kfor... 
    a1for gamma_efor a2for omega_cfor sigma_cfor a3for a4for gcoeffor taxcoeffor... 
    taucfor taulfor taurfor cssfor wssfor labssfor bbarfor ifor_yfor gfor_yfor...
    dss_yfor yssfor mucfor Gwbarfor dss4for Rdssfor Rlssfor omegabarfor Fwbar_wfor... 
    Fwbar_sigfor Gwbar_wfor Gwbar_sigfor sigmawfor ch_y indexss tradess...
    indexforss chfor_y c_y cfor_yfor chss cfssfor cfss chssfor dssEMU yssEMU...
    css cssfor bondss bondforss Pihss Pifss size gamma_ita gamma_itafor... 
    itass itassfor netxss netxforss Fwbar Fwbarfor lpsi lpsifor
elseif mpru==2
    save NI_paramfileopenMPRUFOREIGN betta hcons vartheta theta vareps chi alfa delta Piss Rss Pistar...
    S2p1 Rkss omega_a sigma_a omega_b sigma_b qss rss nss kss bbar_k... 
    a1 gamma_e a2 omega_c sigma_c a3 a4 gcoef taxcoef tauc taul taur css wss labss bbar c_y i_y g_y dss_y yss...
    muc Gwbar gamma_g gamma_c gamma_l gamma_r rhod rhoa rhosigma rhogammae dss4 Rdss Rlss...
    omegabar Fwbar_w Fwbar_sig Gwbar_w Gwbar_sig sigmaw zet varphi varphifor omegaB...
    thetafor Pissfor Rssfor Pistarfor S2p1for Rkssfor omega_afor sigma_afor...
    omega_bfor sigma_bfor qssfor rssfor nssfor kssfor bbar_kfor... 
    a1for gamma_efor a2for omega_cfor sigma_cfor a3for a4for gcoeffor taxcoeffor... 
    taucfor taulfor taurfor cssfor wssfor labssfor bbarfor ifor_yfor gfor_yfor...
    dss_yfor yssfor mucfor Gwbarfor dss4for Rdssfor Rlssfor omegabarfor Fwbar_wfor... 
    Fwbar_sigfor Gwbar_wfor Gwbar_sigfor sigmawfor ch_y indexss tradess...
    indexforss chfor_y c_y cfor_yfor chss cfssfor cfss chssfor dssEMU yssEMU...
    css cssfor bondss bondforss Pihss Pifss size gamma_ita gamma_itafor...
    itass itassfor netxss netxforss Fwbar Fwbarfor lpsi lpsifor
elseif mpru==3
    save NI_paramfileopenMPRUBOTH betta hcons vartheta theta vareps chi alfa delta Piss Rss Pistar...
    S2p1 Rkss omega_a sigma_a omega_b sigma_b qss rss nss kss bbar_k... 
    a1 gamma_e a2 omega_c sigma_c a3 a4 gcoef taxcoef tauc taul taur css wss labss bbar c_y i_y g_y dss_y yss...
    muc Gwbar gamma_g gamma_c gamma_l gamma_r rhod rhoa rhosigma rhogammae dss4 Rdss Rlss...
    omegabar Fwbar_w Fwbar_sig Gwbar_w Gwbar_sig sigmaw zet varphi varphifor omegaB...
    thetafor Pissfor Rssfor Pistarfor S2p1for Rkssfor omega_afor sigma_afor...
    omega_bfor sigma_bfor qssfor rssfor nssfor kssfor bbar_kfor... 
    a1for gamma_efor a2for omega_cfor sigma_cfor a3for a4for gcoeffor taxcoeffor... 
    taucfor taulfor taurfor cssfor wssfor labssfor bbarfor ifor_yfor gfor_yfor...
    dss_yfor yssfor mucfor Gwbarfor dss4for Rdssfor Rlssfor omegabarfor Fwbar_wfor... 
    Fwbar_sigfor Gwbar_wfor Gwbar_sigfor sigmawfor ch_y indexss tradess...
    indexforss chfor_y c_y cfor_yfor chss cfssfor cfss chssfor dssEMU yssEMU...
    css cssfor bondss bondforss Pihss Pifss size gamma_ita gamma_itafor...
    itass itassfor netxss netxforss Fwbar Fwbarfor lpsi lpsifor
elseif mpru==4
save NI_paramfileopenNOMPRU betta hcons vartheta theta vareps chi alfa delta Piss Rss Pistar...
    S2p1 Rkss omega_a sigma_a omega_b sigma_b qss rss nss kss bbar_k... 
    a1 gamma_e a2 omega_c sigma_c a3 a4 gcoef taxcoef tauc taul taur css wss labss bbar c_y i_y g_y dss_y yss...
    muc Gwbar gamma_g gamma_c gamma_l gamma_r rhod rhoa rhosigma rhogammae dss4 Rdss Rlss...
    omegabar Fwbar_w Fwbar_sig Gwbar_w Gwbar_sig sigmaw zet varphi varphifor omegaB...
    thetafor Pissfor Rssfor Pistarfor S2p1for Rkssfor omega_afor sigma_afor...
    omega_bfor sigma_bfor qssfor rssfor nssfor kssfor bbar_kfor... 
    a1for gamma_efor a2for omega_cfor sigma_cfor a3for a4for gcoeffor taxcoeffor... 
    taucfor taulfor taurfor cssfor wssfor labssfor bbarfor ifor_yfor gfor_yfor...
    dss_yfor yssfor mucfor Gwbarfor dss4for Rdssfor Rlssfor omegabarfor Fwbar_wfor... 
    Fwbar_sigfor Gwbar_wfor Gwbar_sigfor sigmawfor ch_y indexss tradess...
    indexforss chfor_y c_y cfor_yfor chss cfssfor cfss chssfor dssEMU yssEMU...
    css cssfor bondss bondforss Pihss Pifss size netxss netxforss...
    gamma_ita gamma_itafor itass itassfor Fwbar Fwbarfor lpsi lpsifor
end