function [ys,params,check] = ERJ_est_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = NK_baseline_steadystate(ys,exo,M_,options_)
% computes the steady state for the NK_baseline.mod and uses a numerical
% solver to do so
% Inputs:
%   - ys        [vector] vector of initial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%   - M_        [structure] Dynare model structure
%   - options   [structure] Dynare options structure
%
% Output:
%   - ys        [vector] vector of steady state values for the the endogenous variables
%   - params    [vector] vector of parameter values
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impose restrictions on parameters)

% Copyright (C) 2013-2020 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
    paramname = M_.param_names{ii};
    eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;

params=M_.params;


%% Enter model equations here

%% The steady state computation

% First initialize known steady state values:
% ss_ERJ_nocbdc;
if phibh == 0
    is_cost=0;
else
    is_cost=1;
end
% keyboard
PSIh=0;
PSIf=0;
PSIw=0;
%import/export shares
omegah_fb=(1-omega_h)*sizef/(sizef+sizew);
omegah_fg=(1-omega_h)*sizef/sizef;
omegah_wb=(1-omega_h)*sizew/(sizef+sizew);
omegaf_hb=(1-omega_f)*sizeh/(sizeh+sizew);
omegaf_hg=(1-omega_f)*sizeh/sizeh;
omegaf_wb=(1-omega_f)*sizew/(sizeh+sizew);
omegaw_hb=(1-omega_w)*sizeh/(sizeh+sizef);
omegaw_fb=(1-omega_w)*sizef/(sizeh+sizef);
omegaw_hg=(1)*sizeh/(sizeh+sizef);
omegaw_fg=(1)*sizef/(sizeh+sizef);
alphah_green=alphah*scale_alpha_green_h;     %technology (green)
alphaf_green=alphaf*scale_alpha_green_f;
% tax_h = 1*scale_taxh + poly;
% tax_f = 1*scale_taxf + poly;
%rhof and rhoh in the steady state
rhoh=rhohss;
rhof=rhofss;
rhow=rhowss;
%Labor in both sectors is calibrated to 1
Lhb=Lsshb;
Lfb=Lssfb;
Lhg=Lsshg;
Lfg=Lssfg;
Lwb=Lsswb;
%MC derives from optimality (i.e. in the equilibrium the mark-up is
%proportional to mrk power) times the tax, nothce there is a different tax
%for each sector

MChb=((nih-1)/(nih))*(1+tauhb);
MCfb=((nif-1)/(nif))*(1+taufb);
MChg=((nih-1)/(nih))*(1+tauhg);
MCfg=((nif-1)/(nif))*(1+taufg);
MCwb=((niw-1)/(niw))*(1+tauwb);
%Exchange rate
NEER=1;
REER=1;
NEERhw=1;
REERhw=1;
NEERfw=1;
REERfw=1;
s=1;
%Risk free rate
rh=1/betah;
rf=1/betaf;
rw=1/betaw;
%All inflation rates are 1 (no inflation)
pih=1;
pif=1;
pihomehb=1;
pihomefb=1;
pihomehg=1;
pihomefg=1;
piforhb=1;
piforfb=1;
piforhg=1;
piforfg=1;
piwrdwb=1;
pihomewb=1;
piforwb=1;
piwrdhb=1;
piwrdfb=1;
pihomewg=1;
piforwg=1;
piw=1;
%New optimal prices (relative to CPI) are also 1 in the ss from F*PO=K and
%F*RER*PO=K
POhomehb=1;
POhomefb=1;
POforhb=1;
POforfb=1;
POhomehg=1;
POhomefg=1;
POforhg=1;
POforfg=1;
POwrdwb=1;
POhomewb=1;
POforwb=1;
POwrdhb=1;
POwrdfb=1;
POhomewg=1;
POforwg=1;

%In the steady state prices are at the optimal level, so...
phomehb=1;
pforfb=1;
pwrdwb=1;
phomehg=1;
pforfg=1;
pforhb=1;
pwrdhb=1;
phomefb=1;
pwrdfb=1;
phomewb=1;
pforwb=1;
pforhg=1;
phomefg=1;
phomewg=1;
pforwg=1;

%Total output brown accounting for the loss in production due to brown
% Ytothb=Khb^(alphah);
% Ytotfb=Kfb^(alphaf);
%This solves the system:
% Y_brown1 = [1-d1*x-d2*x^2]*C1
% C1=A1*L1^(1-alpha)*K1^(alpha)
% x = eta*x(-1) + e1 +e2
% e1 = Y_brown1^(1-gamma)
% e2 = Y_brown2^(1-gamma)
%notice that now also returns on capital and wages are a function of the
%emissions!
% keyboard
% C1=Khb^(alphah);
% C2=Kfb^(alphaf);
%write the function
%give some starting point (no emissions case)
rkhb=((1-betah*(1-deltakh))/(betah));
rkfb=((1-betaf*(1-deltakf))/(betaf));
rkwb=((1-betaw*(1-deltakw))/(betaw));
Khb=(rkhb/(alphah*MChb))^(1/(alphah-1));
Kfb=(rkfb/(alphaf*MCfb))^(1/(alphaf-1));
Kwb=(rkwb/(alphaw*MCwb))^(1/(alphaw-1));
% keyboard
% Eq. 1 Y1 = [1-Loss1]*K1^(alpha1)
% Eq. 2 Y2 = [1-Loss2]*K2^(alpha2)
% Eq. 3 K1 = [(RK1 + tax1*Y1^(1-g1)*(1-g1)*alpha1/K1)/[alpha1*MC1*(1-Loss1)]]^(1/(alpha1-1))
% Eq. 3 K2 = [(RK2 + tax2*Y2^(1-g2)*(1-g2)*alpha2/K2)/[alpha2*MC2*(1-Loss2)]]^(1/(alpha2-1))
% keyboard
% ab1=2;
% ab2=2.8;
% keyboard
F = @(x) [x(1)-(x(3)^alphah)+(dh0+dh1*((1/(1-etaxh))*((1-0)*x(1)^(1-gh)+(1-0)*x(2)^(1-gf)+(1-0)*x(5)^(1-gw)))^dh2)*(x(3)^alphah);
    x(2)-(x(4)^alphaf)+(df0+df1*((1/(1-etaxf))*((1-0)*x(1)^(1-gh)+(1-0)*x(2)^(1-gf)+(1-0)*x(5)^(1-gw)))^df2)*(x(4)^alphaf);
    x(5)-(x(6)^alphaw)+(dw0+dw1*((1/(1-etaxh))*((1-0)*x(1)^(1-gh)+(1-0)*x(2)^(1-gf)+(1-0)*x(5)^(1-gw)))^dw2)*(x(6)^alphaw);
    x(3)-((rkhb+tax_h*((1-0)*x(1)^(1-gh))*(1-gh)*alphah/x(3))/((1-(dh0+dh1*((1/(1-etaxh))*((1-0)*x(1)^(1-gh)+(1-0)*x(2)^(1-gf)+(1-0)*x(5)^(1-gw)))^dh2)) *(alphah*MChb)))^(1/(alphah-1));
    x(4)-((rkfb+tax_f*((1-0)*x(2)^(1-gf))*(1-gf)*alphaf/x(4))/((1-(df0+df1*((1/(1-etaxf))*((1-0)*x(1)^(1-gh)+(1-0)*x(2)^(1-gf)+(1-0)*x(5)^(1-gw)))^df2))*(alphaf*MCfb)))^(1/(alphaf-1));
    x(6)-((rkwb+tax_w*((1-0)*x(5)^(1-gw))*(1-gw)*alphaw/x(6))/((1-(dw0+dw1*((1/(1-etaxh))*((1-0)*x(1)^(1-gh)+(1-0)*x(2)^(1-gf)+(1-0)*x(5)^(1-gw)))^dw2)) *(alphaw*MCwb)))^(1/(alphaw-1))];

% sol=my_system(x,alphah,dh0,dh1,dh2,etaxh,etaxf,gh,gf,alphaf,df0,df1,df2,rkhb,rkfb,tax_h,tax_f,ab1,ab2,MChb,MCfb)
% [x,fval,EXITFLAG] = fsolve(@(x)my_system(x,alphah,dh0,dh1,dh2,etaxh,etaxf,gh,gf,alphaf,df0,df1,df2,rkhb,rkfb,tax_h,tax_f,ab1,ab2,MChb,MCfb) ,[Khb^(alphah),Kfb^(alphaf),Khb,Kfb,0.01,0.01],options);
options = optimoptions('fsolve','Display','none');
%use as initial conditions the solution without emissions effect. The
%solution will be close to that
abh=(tax_h*((Khb^(alphah))^(-gh))/(ab1*ab2))^(1/(ab2-1));
abf=(tax_f*((Kfb^(alphaf))^(-gf))/(ab1*ab2))^(1/(ab2-1));
abw=(tax_w*((Kwb^(alphaw))^(-gw))/(ab1*ab2))^(1/(ab2-1));
try
    [x,fval,EXITFLAG] = fsolve(F,[Khb^(alphah),Kfb^(alphaf),Khb,Kfb,Kwb^(alphaw),Kwb],options);
    if EXITFLAG~=1
        error('Solution for Ytothb and Ytotfb not found');
    end
catch
    check=1;
    return;
end

Ytothb=x(1);
Ytotfb=x(2);
Khb = x(3);
Kfb = x(4);
Ytotwb=x(5);
Kwb = x(6);
abh = 0;
abf = 0;
%compute the steady state environmental cost
eh=(1-abh)*Ytothb^(1-gh);
ef=(1-abf)*Ytotfb^(1-gf);
ew=(1-abw)*Ytotwb^(1-gw);
xh=((1/(1-etaxh)))*(eh+ef+ew);
xf=((1/(1-etaxf)))*(eh+ef+ew);
lossh = dh0 + dh1*xh^dh2;
lossf = df0 + df1*xf^df2;
lossw = dw0 + dw1*xh^dw2;
Lossh = dh0 + dh1*xh^dh2;
Lossf = df0 + df1*xf^df2;
Lossw = dw0 + dw1*xh^dw2;
%Rate of return on risky assets net of loss
rkhb=((1-betah*(1-deltakh))/(betah));
rkfb=((1-betaf*(1-deltakf))/(betaf));
rkwb=((1-betaw*(1-deltakw))/(betaw));
rkhg=((1-betah*(1-deltakh))/(betah));
rkfg=((1-betaf*(1-deltakf))/(betaf));

%compute the to

%Capital (from production function)
%we need to rely again on the non-linear solver
% recall L=1 and Y*sub=t*e, Y=(1-loss)*K^(a)
% rkhg = (1-dh0-dh1*xh^dh2)*(alphah_green*(MChg+sub_h)*exp(Ahg)*(Khg^(alphah_green-1))*Lhg^(1-alphah));
% rkfg = (1-df0-df1*xf^df2)*(alphaf_green*(MCfg+sub_f)*exp(Afg)*(Kfg^(alphaf_green-1))*Lfg^(1-alphaf));
% keyboard
% Eq. 1 Rk = (1-loss1)*alpha_green*MC*K1^(alpha-1) + sub_1*tax_1*e1/(K1)
% rkhg = (1-lossh)*(alphah_green*(MChg+tax_h*eh/((1-lossh)*x(1)^alphah_green))*(Khg^(alphah_green-1)));
% rkfg = (1-lossf)*(alphaf_green*(MCfg+tax_f*ef/((1-lossf)*x(1)^alphaf_green))*(Kfg^(alphaf_green-1)));

%rkhg = (1-lossh)*(alphah_green*(MChg+sub_h)*(Khg^(alphah_green-1)));
% keyboard
F = @(x) [rkhg - (1-lossh)*(alphah_green*(MChg+isrebate*tax_h*eh/((1-lossh)*(x(1)^alphah_green)))*(x(1)^(alphah_green-1)))];
% rkhg = (1-lossh)*(alphah_green*(MChg+sub_h)*(Khg^(alphah_green-1)));
% keyboard

options = optimoptions('fsolve','Display','none','OptimalityTolerance', 1e-10);
%use as initial conditions the solution without emissions effect. The
%solution will be close to that
try
    [x,fval,EXITFLAG] = fsolve(F,[1],options);
    if EXITFLAG~=1
        error('Solution for Ytothb and Ytotfb not found');
    end
catch
    check=1;
    return;
end
Khg=x(1);
F = @(x) [rkfg - (1-lossf)*(alphaf_green*(MCfg+isrebate*tax_f*ef/((1-lossf)*(x(1)^alphaf_green)))*(x(1)^(alphaf_green-1))) ];


options = optimoptions('fsolve','Display','none','OptimalityTolerance', 1e-10);
%use as initial conditions the solution without emissions effect. The
%solution will be close to that
try
    [x,fval,EXITFLAG] = fsolve(F,[Khg],options);
    if EXITFLAG~=1
        error('Solution for Ytothb and Ytotfb not found');
    end
catch
    check=1;
    return;
end

Kfg=x(1);
% Khg=(rkhg/((1-lossh)*alphah_green*MChg))^(1/(alphah_green-1));
% Kfg=(rkfg/((1-lossh)*alphaf_green*MCfg))^(1/(alphaf_green-1));
%Total output green
%Also in this case we have the account for the loss in production due to
%brown
Ytothg=(1-lossh)*Khg^(alphah_green);
Ytotfg=(1-lossf)*Kfg^(alphaf_green);
%compute hte subsidy
sub_h=tax_h*eh*isrebate/Ytothg;
sub_f=tax_f*ef*isrebate/Ytotfg;

%Investments
Ihb=deltakh*Khb;
Ifb=deltakf*Kfb;
Ihg=deltakh*Khg;
Ifg=deltakf*Kfg;
Iwb=deltakw*Kwb;
Iwg=0;
Ih = Ihg + Ihb;
If = Ifg + Ifb;
Iw = Iwb+Iwg;
ratiohb=1;
ratiofb=1;
ratiowb=1;
ratiohg=1;
ratiofg=1;
%Total output of the economy
Ytoth=Ytothb+Ytothg;
Ytotf=Ytotfb+Ytotfg;
Ytotw=Ytotwb;
% brownh=0.9;
% brownf=0.9;

%solve for brown from investments
% I_brown = brown*I_total
% I_green = (1-brown)*I_total
% F = @(x) [Ihb - x(1)*x(2);
%         Ihg - (1-x(1))*x(2)
%           ];
%
% options = optimoptions('fsolve','Display','none','Algorithm','levenberg-marquardt');
% %use as initial conditions the solution without emissions effect. The
% %solution will be close to that
% [x,fval,EXITFLAG] = fsolve(F,[0,0],options);
% F = @(x) [Ifb - x(1)*x(2);
%         Ifg - (1-x(1))*x(2)
%           ];
% brownh=x(1);
% Ih=x(2);
% options = optimoptions('fsolve','Display','none','Algorithm','levenberg-marquardt');
% %use as initial conditions the solution without emissions effect. The
% %solution will be close to that
% [x,fval,EXITFLAG] = fsolve(F,[0,0],options);
% brownf=x(1);
% If=x(2);
%to solve write total output (Ytoth) as Chh+Ghh+Ihh + Exports
%with exports = Chf+Ghf+Ihf+Chw+Ghw+Ihw
%Do this for all countries, you can callect common terms and end up with a
%system of 3 equations and 3 unknowns:
% Ytoth = a1*(Ch+Gh) + a2*(Cf+Gf)+ a3*(Cw+Gw) + Ihh + Ihf + Ihw;
% Ytotf = b1*(Ch+Gh) + b2*(Cf+Gf)+ b3*(Cw+Gw) + Ifh + Iff + Ifw;
% Ytotw = c1*(Ch+Gh) + c2*(Cf+Gf)+ c3*(Cw+Gw) + Iwh + Iwf + Iww;
% weights are simply shares and we can back out investments from total investments
% K1 = Ytoth - (Ihh + Ihf + Ihw);
my_a1=omega_h*brownh+omega_h*(1-brownh);
my_a2=omegaf_hb*brownf+omegaf_hg*(1-brownf);
my_a3=omegaw_hb*brownw+omegaw_hg*(1-brownw);
my_b1=omegah_fb*brownh+omegah_fg*(1-brownh);
my_b2=omega_f*brownf+omega_f*(1-brownf);
my_b3=omegaw_fb*brownw+omegaw_fg*(1-brownw);
my_c1=omegah_wb*brownh;
my_c2=omegaf_wb*brownf;
my_c3=omega_w*brownw;
my_K1=Ytoth - (omega_h*Ihb+omega_h*Ihg+omegaf_hb*Ifb+omegaw_hb*Iwb+omegaf_hg*Ifg);
my_K2=Ytotf - (omega_f*Ifb+omega_f*Ifg+omegah_fb*Ihb+omegaw_fb*Iwb+omegah_fg*Ihg);
my_K3=Ytotw - (omega_w*Iwb+omegah_wb*Ihb+omegaf_wb*Ifb);
%auxiliary variables
my_psi=my_K2-(my_b1/my_a1)*my_K1;
my_alpha=my_b2-my_b1*my_a2/my_a1;
my_beta=my_b3-my_b1*my_a3/my_a1;
tmp1 = my_K3 - my_c1*my_K1/my_a1 + my_c1*my_a2*my_psi/(my_a1*my_alpha) - my_c2*my_psi/my_alpha;
tmp2 = my_c1*my_a2*my_beta/(my_a1*my_alpha) - my_c1*my_a3/my_a1 - my_c2*my_beta/my_alpha + my_c3;
CGw = tmp1/tmp2;
CGf = (my_psi-my_beta*CGw)/my_alpha;
CGh = (my_K1 - my_a2*CGf - my_a3*CGw)/my_a1;
%alternative solution with symbolic solver
% syms x1 x2 x3
% eqns = [my_K1 == my_a1*x1+my_a2*x2+my_a3*x3 , my_K2 == my_b1*x1+my_b2*x2+my_b3*x3, my_K3 == my_c1*x1+my_c2*x2+my_c3*x3];
% S = solve(eqns,[x1 x2 x3]);

%now get public spending
%Define G from the share of total output
Gh=GYhss*Ytoth;
Gf=GYfss*Ytotf;
Gw=GYwss*Ytotw;
%Define all subcomponents of G
%Domestic brown
Gh_hb = omega_h*brownh*Gh;
Gf_fb = omega_f*brownf*Gf;
Gw_wb = omega_w*brownw*Gw;
%Domestic green
Gh_hg = omega_h*(1-brownh)*Gh;
Gf_fg = omega_f*(1-brownf)*Gf;
%Imported brown
Gh_fb = (omegah_fb)*brownh*Gh;
Gh_wb = (omegah_wb)*brownh*Gh;
Gf_hb = (omegaf_hb)*brownf*Gf;
Gf_wb = (omegaf_wb)*brownf*Gf;
Gw_hb = (omegaw_hb)*brownw*Gw;
Gw_fb = (omegaw_fb)*brownw*Gw;

%Imported green
Gh_fg = (omegah_fg)*(1-brownh)*Gh;
Gf_hg = (omegaf_hg)*(1-brownf)*Gf;
Gw_hg = (omegaw_hg)*(1-brownw)*Gw;
Gw_fg = (omegaw_fg)*(1-brownw)*Gw;
Ih=Ihb+Ihg;
If=Ifb+Ifg;
Iw=Iwb+Iwg;
%once we have G we can get C
Ch = CGh - Gh;
Cf = CGf - Gf;
Cw = CGw - Gw;

%demand for domestic goods (brown) (52) (53) (54)
Ch_hb = omega_h*brownh*Ch;
Cf_fb = omega_f*brownf*Cf;
Cw_wb = omega_w*brownw*Cw;

%demand for domestic goods (green) (55) (56)
Ch_hg = omega_h*(1-brownh)*Ch;
Cf_fg = omega_f*(1-brownf)*Cf;

%demand for imported goods (brown) (57) (58) (59) (60) (61) (62)
Ch_fb = (omegah_fb)*brownh*Ch;
Ch_wb = (omegah_wb)*brownh*Ch;
Cf_hb = (omegaf_hb)*brownf*Cf;
Cf_wb = (omegaf_wb)*brownf*Cf;
Cw_hb = (omegaw_hb)*brownw*Cw;
Cw_fb = (omegaw_fb)*brownw*Cw;

%demand for imported goods (green) (63) (64) (65) (66)
Ch_fg = (omegah_fg)*(1-brownh)*Ch;
Cf_hg = (omegaf_hg)*(1-brownf)*Cf;
Cw_hg = (omegaw_hg)*(1-brownw)*Cw;
Cw_fg = (omegaw_fg)*(1-brownw)*Cw;

%fix investments
Ih=Ihb+Ihg;
If=Ifb+Ifg;
Iw=Iwb+Iwg;
% % compute brown share based on Investments
Ih_hb = omega_h*Ihb;
If_fb = omega_f*Ifb;
Iw_wb = omega_w*Iwb;

%Consumption of domestic investment goods (green) (88) (89)
Ih_hg = omega_h*Ihg;
If_fg = omega_f*Ifg;

%consumption of imported investment goods (brown) (90) (91) (92) (93) (94) (95)
Ih_fb = (omegah_fb)*Ihb;
Ih_wb = (omegah_wb)*Ihb;

If_hb = (omegaf_hb)*Ifb;
If_wb = (omegaf_wb)*Ifb;

Iw_hb = (omegaw_hb)*1*Iwb;
Iw_fb = (omegaw_fb)*1*Iwb;

%consumption of imported investment goods (green) (96) (97) (98) (99)
Ih_fg = (omegah_fg)*Ihg;
If_hg = (omegaf_hg)*Ifg;
Iw_hg=0;
Iw_fg=0;
%now re-aggregate into domestic and exports
Ybh = Ch_hb + Gh_hb + Ih_hb;
Ybf = Cf_fb + Gf_fb + If_fb;
Ybw = Cw_wb + Gw_wb + Iw_wb;
Ygh = Ch_hg + Gh_hg + Ih_hg;
Ygf = Cf_fg + Gf_fg + If_fg;
Yh_forb = Cf_hb + Gf_hb + If_hb;
Yh_forbw = Cw_hb + Gw_hb + Iw_hb;
Yf_forb = Ch_fb + Gh_fb + Ih_fb;
Yf_forbw = Cw_fb + Gw_fb + Iw_fb;
Yw_forbh = Ch_wb + Gh_wb + Ih_wb;
Yw_forbf = Cf_wb + Gf_wb + If_wb;
Yh_forg = Cf_hg + Gf_hg + If_hg;
Yh_forgw = Cw_hg + Gw_hg;
Yf_forg = Ch_fg + Gh_fg + Ih_fg;
Yf_forgw = Cw_fg + Gw_fg;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Old calculatins starting here
% CGh = Ytoth - Ihb -Ihg;
% CGf = Ytotf - Ifb -Ifg;
% CGw = Ytotw - Iwb -Iwg;
%
% Ybh =omega_h*brownh*CGh + omega_h*Ihb;
% Ybf =omega_f*brownf*CGf + omega_f*Ifb;
% Ybw =omega_w*brownw*CGw + omega_w*Iwb;
%
% Ygh =omega_h*(1-brownh)*CGh+omega_h*Ihg;
% Ygf =omega_f*(1-brownf)*CGf+omega_f*Ifg;
% Ygw =0;
%
% %exports of h country (brown)
% Yh_forb =(omegaf_hb)*brownf*CGf+(omegaf_hb)*Ifb;
% Yh_forbw=(omegaw_hb)*brownw*CGw+(omegaw_hb)*Iwb;
%
% %exports of f country (brown)
% Yf_forb =(omegah_fb)*brownh*CGh+(omegah_fb)*Ihb;
% Yf_forbw=(omegaw_fb)*brownw*CGw+(omegaw_fb)*Iwb;
%
% %export of w country (brown)
% Yw_forbh=(omegah_wb)*brownh*CGh+(omegah_wb)*Ihb;
% Yw_forbf=(omegaf_wb)*brownf*CGf+(omegaf_wb)*Ifb;
%
% %exports of h country (green)
% Yh_forg =(omegaf_hg)*(1-brownf)*CGf+(omegaf_hg)*Ifg;
% Yh_forgw=(omegaw_hg)*(1-brownw)*CGw+(omegaw_hg)*Iwg;
%
% %exports of h country (green)
% Yf_forg =(omegah_fg)*(1-brownh)*CGh+(omegah_fg)*Ihg;
% Yf_forgw=(omegaw_fg)*(1-brownw)*CGw+(omegaw_fg)*Iwg;
% % Ybh =omegah*brownh*Ytoth;
% % Ybf =omegaf*brownf*Ytotf;
% % Ygh =omegah*(1-brownh)*Ytoth;
% % Ygf =omegaf*(1-brownf)*Ytotf;
% % Yh_forb =(1-omegaf)*brownf*Ytotf;
% % Yf_forb =(1-omegah)*brownh*Ytoth;
% % Yh_forg =(1-omegaf)*(1-brownf)*Ytotf;
% % Yf_forg =(1-omegah)*(1-brownh)*Ytoth;
%
%
% %Define G from the share of total output
% Gh=GYhss*Ytoth;
% Gf=GYfss*Ytotf;
% Gw=GYwss*Ytotw;
% %Define all subcomponents of G
% %Domestic brown
% Gh_hb = omega_h*brownh*Gh;
% Gf_fb = omega_f*brownf*Gf;
% Gw_wb = omega_w*brownw*Gw;
% %Domestic green
% Gh_hg = omega_h*(1-brownh)*Gh;
% Gf_fg = omega_f*(1-brownf)*Gf;
% %Imported brown
% Gh_fb = (omegah_fb)*brownh*Gh;
% Gh_wb = (omegah_wb)*brownh*Gh;
% Gf_hb = (omegaf_hb)*brownf*Gf;
% Gf_wb = (omegaf_wb)*brownf*Gf;
% Gw_hb = (omegaw_hb)*brownw*Gw;
% Gw_fb = (omegaw_fb)*brownw*Gw;
%
% %Imported green
% Gh_fg = (omegah_fg)*(1-brownh)*Gh;
% Gf_hg = (omegaf_hg)*(1-brownf)*Gf;
% Gw_hg = (omegaw_hg)*(1-brownw)*Gw;
% Gw_fg = (omegaw_fg)*(1-brownw)*Gw;
% Ih=Ihb+Ihg;
% If=Ifb+Ifg;
% Iw=Iwb+Iwg;
% %Define all subcomponents of I
% % compute brown share based on Investments
% % Ih_hb = omegah*brownh*Ih;
% % If_fb = omegaf*brownf*If;
% % Ih_hg = omegah*(1-brownh)*Ih;
% % If_fg = omegaf*(1-brownf)*If;
% % %Import brown
% % Ih_fb = (1-omegah)*brownh*Ih;
% % If_hb = (1-omegaf)*brownf*If;
% % %Import green
% % Ih_fg = (1-omegah)*(1-brownh)*Ih;
% % If_hg = (1-omegaf)*(1-brownf)*If;
%
% % % compute brown share based on Investments
% Ih_hb = omega_h*Ihb;
% If_fb = omega_f*Ifb;
% Iw_wb = omega_w*Iwb;
%
% %Consumption of domestic investment goods (green) (88) (89)
% Ih_hg = omega_h*Ihg;
% If_fg = omega_f*Ifg;
%
% %consumption of imported investment goods (brown) (90) (91) (92) (93) (94) (95)
% Ih_fb = (omegah_fb)*Ihb;
% Ih_wb = (omegah_wb)*Ihb;
%
% If_hb = (omegaf_hb)*Ifb;
% If_wb = (omegaf_wb)*Ifb;
%
% Iw_hb = (omegaw_hb)*1*Iwb;
% Iw_fb = (omegaw_fb)*1*Iwb;
%
% %consumption of imported investment goods (green) (96) (97) (98) (99)
% Ih_fg = (omegah_fg)*Ihg;
% If_hg = (omegaf_hg)*Ifg;
% Iw_hg=0;
% Iw_fg=0;

% %Now turn to consumption. Total consumption is
% % Ybh = omegah*brownh*C+omegah*brownh*I+omegah*brownh*G
% % Ybf = omegaf*brownf*C+omegaf*brownf*I+omegaf*brownf*G
% % Ch = (Ybh - omegah*brownh*Ih - omegah*brownh*Gh)/(omegah*brownh);
% % Cf = (Ybf - omegaf*brownh*If - omegaf*brownf*Gf)/(omegaf*brownf);
% Ch=CGh-Gh;
% Cf=CGf-Gf;
% Cw=CGw-Gw;
% %recover sub components of consumption
% %demand for domestic goods (brown) (52) (53) (54)
% Ch_hb = omega_h*brownh*Ch;
% Cf_fb = omega_f*brownf*Cf;
% Cw_wb = omega_w*brownw*Cw;
%
% %demand for domestic goods (green) (55) (56)
% Ch_hg = omega_h*(1-brownh)*Ch;
% Cf_fg = omega_f*(1-brownf)*Cf;
%
% %demand for imported goods (brown) (57) (58) (59) (60) (61) (62)
% Ch_fb = (omegah_fb)*brownh*Ch;
% Ch_wb = (omegah_wb)*brownh*Ch;
% Cf_hb = (omegaf_hb)*brownf*Cf;
% Cf_wb = (omegaf_wb)*brownf*Cf;
% Cw_hb = (omegaw_hb)*brownw*Cw;
% Cw_fb = (omegaw_fb)*brownw*Cw;
%
% %demand for imported goods (green) (63) (64) (65) (66)
% Ch_fg = (omegah_fg)*(1-brownh)*Ch;
% Cf_hg = (omegaf_hg)*(1-brownf)*Cf;
% Cw_hg = (omegaw_hg)*(1-brownw)*Cw;
% Cw_fg = (omegaw_fg)*(1-brownw)*Cw;
%
%
% % Ybh = (Ch_hb + Gh_hb + Ih_hb)
% % Ybf =Cf_fb + Gf_fb + If_fb;
% % Ygh =Ch_hg + Gh_hg + Ih_hg;
% % Ygf =Cf_fg + Gf_fg + If_fg;
% % Yh_forb =Cf_hb + Gf_hb + If_hb;
% % Yf_forb =Ch_fb + Gh_fb + Ih_fb;
% % Yh_forg =Cf_hg + Gf_hg + If_hg;
% % Yf_forg =Ch_fg + Gh_fg + Ih_fg;
% % Ytoth = Ygh + Yh_forg + Ybh + Yh_forb
% %
% % Ytothg = Ygh+Yh_forg + agg1g;
% % Ytotfg = Ygf+Yf_forg + agg2g;
% % Ytothb = Ybh+Yh_forb + agg1b;
% % Ytotfb = Ybf+Yf_forb + agg2b;
%Old calculations ending here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%real wages brown
% whb = (1-dh0-dh1*xh^dh2)*(1-alphah)*MChb*exp(Ahb)*(Khb^(alphah))*Lhb^(-alphah) - tax_h*(1-gh)*eh*(1-alphah)/Lhb
whb=(1-lossh)*(1-alphah)*MChb*Khb^(alphah)-tax_h*eh*(1-gh)*(1-alphah);
wfb=(1-lossf)*(1-alphaf)*MCfb*Kfb^(alphaf)-tax_f*ef*(1-gf)*(1-alphaf);
wwb=(1-lossw)*(1-alphaw)*MCwb*((Kwb)^(alphaw))*Lwb^(-alphaw) - tax_w*(1-gw)*ew*(1-alphaw)/Lwb;

%real wages green
whg=(1-lossh)*(1-alphah)*(MChg+sub_h)*Khg^(alphah_green);
wfg=(1-lossf)*(1-alphaf)*(MCfg+sub_f)*Kfg^(alphaf_green);

%get lambda from Euler equation
Lambdah =((Ch - hh*Ch)^(-sigmah)) - betah*hh*((Ch - hh*Ch)^(-sigmah));
Lambdaf =((Cf - hf*Cf)^(-sigmaf)) - betaf*hf*((Cf - hf*Cf)^(-sigmaf));
Lambdaw =((Cw - hw*Cw)^(-sigmaw)) - betaw*hw*((Cw - hw*Cw)^(-sigmaw));

%Definition of Q
Qhb=Lambdah;
Qfb=Lambdaf;
Qhg=Lambdah;
Qfg=Lambdaf;
Qwb=Lambdaw;

Ghss=GYhss*Ytoth;
Gfss=GYfss*Ytotf;
Gwss=GYwss*Ytotw;

%Items of price setting (brown)
FFhb=Lambdah*Ybh/(1-betah*csih);
KKhb=FFhb;
FFfb=Lambdaf*Ybf/(1-betaf*csif);
KKfb=FFfb;
FFwb=Lambdaw*Ybw/(1-betaw*csiw);
KKwb=FFwb;

FFh_forb = Lambdah*Yh_forb/(1-betah*csih);
FFh_forwb=Lambdah*Yh_forbw/(1-betah*csih);

FFf_forb = Lambdaf*Yf_forb/(1-betaf*csif);
FFf_forwb= Lambdaf*Yf_forbw/(1-betaf*csif);

FFf_forwhb = Lambdaw*Yw_forbh/(1-betaw*csiw);
FFf_forwfb = Lambdaw*Yw_forbf/(1-betaw*csiw);


KKh_forb=FFh_forb;
KKf_forb=FFf_forb;
KKh_forwb=FFh_forwb;
KKf_forwb=FFf_forwb;
KKf_forwhb=FFf_forwhb;
KKf_forwfb=FFf_forwfb;
%Items of price setting (green)
FFhg=Lambdah*Ygh/(1-betah*csih);
KKhg=FFhg;
FFfg=Lambdaf*Ygf/(1-betaf*csif);
KKfg=FFfg;
FFh_forg = Lambdah*Yh_forg/(1-betah*csih);
FFf_forg = Lambdaf*Yf_forg/(1-betaf*csif);
KKh_forg = FFh_forg;
KKf_forg = FFf_forg;
FFh_forwg = Lambdah*Yh_forgw/(1-betah*csih);
FFf_forwg = Lambdaf*Yf_forgw/(1-betaf*csif);
KKh_forwg=FFh_forwg;
KKf_forwg=FFf_forwg;

%forign bond holdings
if is_cost
    Bh_for=(betah*rf-1)/phibh;
    Bf_for=(betaf*rh-1)/phibf;
    Bh_for_w=(betah*rw-1)/phibh;
    Bf_for_w=(betaf*rw-1)/phibf;
    Bw_for_h=(betaw*rh-1)/phibw;
    Bw_for_f=(betaw*rf-1)/phibw;
else
    Bh_for=0;
    Bf_for=0;
    Bh_for_w=0;
    Bf_for_w=0;
    Bw_for_h=0;
    Bw_for_f=0;
end
% Bh_for=(-Lambdah + rf)/(Lambdah*phibh);
% Bf_for= (-Lambdaf  + rh )/(Lambdaf*phibf);

%domestic bond holdings
Bh=-(sizef*Bf_for+sizew*Bw_for_h)/sizeh;
Bf=-(sizeh*Bh_for+sizew*Bw_for_f)/sizef;
Bw=-(sizeh*Bh_for_w + sizef*Bf_for_w)/sizew;
% Bh=-Bf_for ;
% Bf=-Bh_for;

%price dispersion (all 1 in ss)
dhomehb = 1;
dforhb  = 1;
dforhbw = 1;
dhomefb = 1;
dforfb  = 1;
dforfbw = 1;
dwb     = 1;
dwb2    = 1;
dwb3    = 1;
dhomehg = 1;
dhomefg = 1;
dforhgw = 1;
dforhg  = 1;
dforfg  = 1;
dforfgw = 1;


%Steady state tax
% Tss=Bh + NEER*Bh_for + Ch + Ihb + Ihg - rh*Bh- NEER*rf*Bh_for +(phibh/2)*((NEER*Bh_for)^2) - rkhb*Khb - rkhg*Khg - whb*Lhb - whg*Lhg-((phomehg - MChg*dhomehg)*Ygh + (phomefg*REER - MChg*dforhg)*Yh_forg + (phomehb - MChb*dhomehb)*Ybh + (phomefb*REER - MChg*dforhb)*Yh_forb);
%Tss2 = (Bh+REER*Bh_for+Ch+Ihb+Ihg) - (rh*Bh/pih+REER*rf*Bh_for/pif-(phibh/2)*((REER*Bh_for)^2)+rkhb*Khb+rkhg*Khg+whb*Lhb+whg*Lhg+(phomehg - MChg*dhomehg)*Ygh + (phomefg*REER - MChg*dforhg)*Yh_forg + (phomehb - MChb*dhomehb)*Ybh + (phomefb*REER - MChg*dforhb)*Yh_forb);
tmp1=(Bh+REER*Bh_for+REERhw*Bh_for_w+Ch+Ihb+Ihg);
tmp2=rh*Bh + REER*rf*Bh_for + REERhw*rw*Bh_for_w-(phibh/2)*((REER*Bh_for)^2) - (phibh/2)*((REERhw*Bh_for_w)^2) + ...
    rkhb*Khb + rkhg*Khg + whb*Lhb + whg*Lhg  + ...
    (sizeh/sizeh)*(phomehg - MChg*dhomehg)*Ygh + (sizef/sizeh)*(phomefg*REER - MChg*dforhg)*Yh_forg + ...
    (sizeh/sizeh)*(phomehb - MChb*dhomehb)*Ybh + (sizef/sizeh)*(phomefb*REER - MChb*dforhb)*Yh_forb + ...
    (sizew/sizeh)*(phomewb*REERhw - MChb*dforhbw)*Yh_forbw + (sizew/sizeh)*(phomewg*REERhw - MChg*dforhgw)*Yh_forgw;
Tss = tmp1-tmp2;

tmp1=(Bf+Bf_for/REER+REERfw*Bf_for_w+Cf+Ifb+Ifg);
tmp2=rf*Bf + (1/REER)*rh*Bf_for + REERfw*rw*Bf_for_w- (phibf/2)*((Bf_for/REER)^2) - (phibf/2)*((REERfw*Bf_for_w)^2) + ...
    rkfb*Kfb + rkfg*Kfg + wfb*Lfb + wfg*Lfg + ...
    (sizef/sizef)*(pforfg - MCfg*dhomefg)*Ygf + (sizeh/sizef)*(pforhg/REER - MCfg*dforfg)*Yf_forg + ...
    (sizef/sizef)*(pforfb - MCfb*dhomefb)*Ybf + (sizeh/sizef)*(pforhb/REER - MCfb*dforfb)*Yf_forb + ...
    (sizew/sizef)*(pforwb*REERfw - MCfb*dforfbw)*Yf_forbw + (sizew/sizef)*(pforwg*REERfw - MCfg*dforfgw)*Yf_forgw;
Tfss=tmp1-tmp2;

%shocks
ECh = 0;
ECf = 0;
ECw = 0;
Ahb = 0;
Afb = 0;
Awb = 0;
RRh = 0;
RRf = 0;
RRw = 0;
Ph  = 0;
Pf  = 0;
Pw  = 0;

%labor is a bit tricky as one needs to define 3 variables chihb chihg Lh
%using these 3 equations
%  chihb * (Lh^(phih+sigmalh))*Lhb^(sigmalh) = whb*Lambdah;
%  chihg * (Lh^(phih+sigmalh))*Lhg^(sigmalh) = whg*Lambdah;
%  Lh=(chihb*Lhb^(1+sigmalh) + chihg*Lhg^(1+sigmalh))^(1/(1+sigmalh));
%  solve the first two for chihb and chihg
%  chihb = (Lh^(-phih-sigmalh))*whb*Lambdah / (Lhb^(sigmalh))
%  chihg = (Lh^(-phih-sigmalh))*whg*Lambdah / (Lhg^(sigmalh))
%  call tmp1 == whb*Lambdah / (Lhb^(sigmalh)) ; tmp2 == whg*Lambdah / (Lhg^(sigmalh))
%  call tmp3 == Lhb^(1+sigmalh) ; tmp4 == Lhg^(1+sigmalh)
%  rewrite the third equation
%  Lh^(1+sigmalh) = chihb*tmp3 + chihg*tmp4
%  substitute
%  Lh^(1+sigmalh) = tmp1*tmp3*(Lh^(-phih-sigmalh)) + tmp2*tmp4*(Lh^(-phih-sigmalh))
%  collect Lh
%  Lh^(1+sigmalh+phih+sigmalh) = tmp1*tmp3 + tmp2*tmp4;
%  Lh= [tmp1*tmp3 + tmp2*tmp4]^(1/(1+sigmalh+phih+sigmalh))
tmp1 = whb*Lambdah / (Lhb^(sigmalh));
tmp2 = whg*Lambdah / (Lhg^(sigmalh));
tmp3 = Lhb^(1+sigmalh);
tmp4 = Lhg^(1+sigmalh);
Lh= [tmp1*tmp3 + tmp2*tmp4]^(1/(1+sigmalh+phih-sigmalh));
chihb = whb*Lambdah/((Lh^(phih-sigmalh))*Lhb^(sigmalh));
chihg = whg*Lambdah/((Lh^(phih-sigmalh))*Lhg^(sigmalh));

tmp1 = wfb*Lambdaf / (Lfb^(sigmalf));
tmp2 = wfg*Lambdaf / (Lfg^(sigmalf));
tmp3 = Lfb^(1+sigmalf);
tmp4 = Lfg^(1+sigmalf);
Lf= [tmp1*tmp3 + tmp2*tmp4]^(1/(1+sigmalf+phif-sigmalf));
chifb = wfb*Lambdaf/((Lf^(phif-sigmalf))*Lfb^(sigmalf));
chifg = wfg*Lambdaf/((Lf^(phif-sigmalf))*Lfg^(sigmalf));

Lwg=0;

chiwb = (wwb*Lambdaw)^((1+sigmalw)/(1+phiw));
chiwg=1;
Lw = (chiwb*Lwb^(1+sigmalw) + chiwg*Lwg^(1+sigmalw))^(1/(1+sigmalw));
% if 0
% %This only if calvo pricing on wages
% chih=Lambdah*((niwh-1)/niwh)*wh;
% chif=Lambdaf*((niwf-1)/niwf)*wf;
%
% FWh = Lh*Lambdah*((niwh-1)/niwh)/(1-betah*csiwh);
% FWf = Lf*Lambdaf*((niwf-1)/niwf)/(1-betaf*csiwf);
% KWh = (chih*(Lh^(1+phih))*wh^(phih*niwh))/(1-betah*csiwh);
% KWf = (chif*(Lf^(1+phif))*wf^(phif*niwf))/(1-betaf*csiwf);
%
% woh = (KWh/FWh)^(1/(1+niwh*phih));
% wof = (KWf/FWf)^(1/(1+niwf*phif));
% end

% %Utility
% Uh = ((exp(ECh)*log(Ch - hh*Ch)) ) - ((1/(1+phih))*Lh^(1+phih))  ;
% Uf = ((exp(ECf)*log(Cf - hf*Cf)) ) - ((1/(1+phif))*Lf^(1+phif))  ;
%
% %Welfare
% Welh=Uh/(1-betah);
% Welf=Uf/(1-betaf);
% Wh=Welh;
% Wf=Welf;
% % ForDemh=ForDemh_ss;
% % ForDemf=ForDemf_ss;
% Ph=0;
% Pf=0;
%Accounting for numerical discrepancies
% agg1b = Ytothb - Ybh - Yh_forb;
% agg2b = Ytotfb - Ybf - Yf_forb;
% agg1g = Ytothg - Ygh - Yh_forg;
% agg2g = Ytotfg - Ygf - Yf_forg;
% agg1b = (Ytoth - ((sizeh/sizeh)*dhomehb*Ybh + (sizef/sizeh)*dforhb*Yh_forb + (sizew/sizeh)*dforhbw*Yh_forbw +...
%         (sizeh/sizeh)*dhomehg*Ygh*(1) + (sizef/sizeh)*dforhg*Yh_forg*(1) + (sizew/sizeh)*dforhgw*Yh_forgw*(1))/2;
%         agg1b + agg1g;
% Ytotf = (sizef/sizef)*dhomefb*Ybf + (sizeh/sizef)*dforfb*Yf_forb + (sizew/sizef)*dforfbw*Yf_forbw+
%         (sizef/sizef)*dhomefg*Ygf + (sizeh/sizef)*dforfg*Yf_forg + (sizew/sizef)*dforfgw*Yf_forgw+
%         agg2b + agg2g;
Ytoth = dhomehb*Ybh + dforhb*Yh_forb + dforhbw*Yh_forbw +...
    dhomehg*Ygh + dforhg*Yh_forg + dforhgw*Yh_forgw;
Ytotf = dhomefb*Ybf + dforfb*Yf_forb + dforfbw*Yf_forbw+...
    dhomefg*Ygf + dforfg*Yf_forg + dforfgw*Yf_forgw;
Ytotw = dwb*Ybw + dwb2*Yw_forbh + dwb3*Yw_forbf;

agg1b = 0;
agg2b = 0;
agg1g = 0;
agg2g = 0;
%utility
Uh = ((exp(ECh)*log(Ch - hh*Ch)) ) - ((1/(1+phih))*Lh^(1+phih)) ;
Uf = ((exp(ECf)*log(Cf - hf*Cf)) ) - ((1/(1+phif))*Lf^(1+phif)) ;
Uw = ((exp(ECw)*log(Cw - hw*Cw)) ) - ((1/(1+phiw))*Lw^(1+phiw)) ;
%welfare
Wh= Uh/(1-betah);
Wf= Uf/(1-betaf);
Ww= Uw/(1-betaw);
yh_o=0;
yf_o=0;
ch_o=0;
cf_o=0;
ih_o=0;
if_o=0;
pih_o=0;
pif_o=0;
lh_o=0;
lf_o=0;
rh_o=log(rh);
rf_o=log(rf);
fx_o=0;
mf_o=0;
mh_o=0;
lyh=log(Ytoth);
lyf=log(Ytotf);
lneer=log(NEER);
tax_imp_brown_h=0;
tax_imp_brown_f=0;
yh_o=0;
yf_o=0;
pih_o = 0;
pif_o = 0;
rh_o=400*log(rh);
rf_o=400*log(rf);
x_hh = 0;
x_ff = 0;

yh_o=0;
yf_o=0;
ch_o=0;
cf_o=0;
ih_o=0;
if_o=0;
pih_o=0;
pif_o=0;
lh_o=0;
lf_o=0;
rh_o=log(rh);
rf_o=log(rf);
fx_o=0;
mf_o=0;
mh_o=0;
lyh=log(Ytoth);
lyf=log(Ytotf);
lneer=log(NEER);
tax_imp_brown_h=0;
tax_imp_brown_f=0;
yh_o=0;
yf_o=0;
pih_o = 0;
pif_o = 0;
rh_o=400*log(rh);
rf_o=400*log(rf);
x_hh = 0;
x_ff = 0;

Dom_h=(Ybh+Ygh);
Exp_h=(Yh_forb+Yh_forg);
Dom_f=(Ybf+Ygf);
Exp_f=(Yf_forb+Yf_forg);
WWh=whg*(Lhg/(Lhg+Lhb))+whb*(Lhb/(Lhg+Lhb));
Kh=Khb+Khg;
Kf=Khb+Khg;
imp_tax_h =1;
imp_tax_f=1;
imp_tax_hw=1;
imp_tax_fw=1;
imp_tax_wh=1;
imp_tax_wf=1;
imp_tax_w=1;
WW=sizeh*Wh+sizef*Wf+sizew*Uw;
UU=sizeh*Uh+sizef*Uf+sizew*Uw;
PSIhp1=PSIh;
PSIfp1=PSIf;
Ihgm1=Ihg;
Ifgm1=Ifg;
Ihbm1=Ihb;
Ifbm1=Ifb;
Wtot=0.5*Wh+0.5*Wf;
yh_o=100*(log(Ytoth)-log(Ytoth));
yf_o=100*(log(Ytotf)-log(Ytotf));
ch_o=100*(log(Ch)-log(Ch));
cf_o=100*(log(Cf)-log(Cf));
pih_o = 100*(pih-1);
pif_o = 100*(pif-1);
lh_o=100*(log(Lh)-log(Lh));
lf_o=100*(log(Lf)-log(Lf));
rh_o=100*(rh-1/betah);
rf_o=100*(rf-1/betaf);
x_hh = (log(eh)-log(eh));
x_ff = (log(ef)-log(ef));
ih_o=100*(log(Ih)-log(Ih));
if_o=100*(log(If)-log(If));
gh_o=100*(log(Gh)-log(Gh));
gf_o=100*(log(Gf)-log(Gf));
mm_h=rh_o;
mm_f=rf_o;
dmm_h=0;
dmm_f=0;
rh2_o=0;
rf2_o=0;
fx=0;
tot=1;
rh_o2=0;
rf_o2=0;
%% end own model equations

params=NaN(NumberOfParameters,1);
for iter = 1:length(M_.params) %update parameters set in the file
    eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
    varname = M_.endo_names{ii};
    eval(['ys(' int2str(ii) ') = ' varname ';']);
end
