% BS paper: model non-linearised, letting dynare do it
% trying to do it with log-linearisation

@#define countries = [ "1", "2" ]
% Period of news
@#define q = 0
q= @{q}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% Declaring variables  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

var p rerCONS p2T c_rel y_rel;
@#for co in countries
var c@{co} n@{co} k@{co} y@{co} i@{co} zT@{co} zN@{co} qa@{co} qb@{co} a@{co} b@{co} mi@{co} lamda@{co} phi@{co}
cT@{co} cN@{co} yT@{co} yT_I@{co} yN@{co} pN@{co} S@{co} ST@{co} SN@{co} nT@{co} nN@{co} P@{co} u@{co} y@{co}_real
% all new LMs
ksi@{co} pi@{co} zita@{co} tau@{co} deltaU@{co} deltaprimeU@{co} sI@{co} sprimeI@{co}; %nx1_y1 

varexo eT@{co} eN@{co} ;
@#endfor
var rerNT rerLM cN1_cT1 n_rel mi_phi ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% Declaring Parameters  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

@#for co in countries
parameters 
alphaT_@{co} alphaN_@{co} beta_@{co} delta_@{co} gamma_@{co} theta_@{co} rho_@{co} omega_@{co} 
omegaT_@{co} is_@{co} yT_y@{co} yN_y@{co} y@{co}_S@{co} yT@{co}_ST@{co} yN@{co}_SN@{co} c@{co}_y@{co} i@{co}_y@{co}  n_ss_@{co} 
n@{co}_y@{co} ni_@{co} PSI_@{co} kappa_@{co} mi@{co}_phi@{co} qa@{co}_ss qb@{co}_ss 
A@{co}T_11 A@{co}T_12 A@{co}T_13 A@{co}T_14 A@{co}N_11 A@{co}N_12 A@{co}N_13 A@{co}N_14 
sigmaT_@{co} sigmaN_@{co} nN@{co}_yN@{co} pNyN@{co}_SN@{co} 
nT@{co}_yT@{co} p@{co}N_ss cN_c@{co} cT_c@{co} cT_y@{co} ksi@{co}_pi@{co}
P@{co}_ss P@{co}_ss1 u_ss_@{co} elasticity_deltaprime_@{co} deltadoubleprime_@{co} 
deltaprime_@{co} epsilon_@{co} ;
@#endfor
parameters 
p2T_ss
correT % tradeables, across countries
correN  % non-tradeables, across countries
correTN1 % T and NT, country one
correTN2 % T and NT, country two
correTNa % T and NT, across countries --> corr(eT1,eN2)
correTNb ; % T and NT, across countries --> corr(eT2, eN1)

@#for co in countries
% Calibrates values - those that I FIXED to be so. Based on those, I estimate the ss.
alphaT_@{co} =  (1-0.61) ; % Following Corsetti et al
alphaN_@{co} =  (1-0.56) ; % Following Corsetti et al
%alphaT_@{co} =  (1-0.67) ; % Following BT. They do not distinguish between sectors, and set that to the standard RBC calibrated value. Note that by doing so, they estimate Solow residuals accordingly. 
%alphaN_@{co} =  alphaT_@{co} ;
beta_@{co} = 1/1.04; % calibration for annual data, Beningo and Thoenissen (JR is 0.985 for q-data)
delta_@{co} = 0.1;
gamma_@{co} = -1;
theta_@{co} = 2; 
epsilon_@{co} = 2;
rho_@{co} = 0.74;
ni_@{co} = 1.4 ;
kappa_@{co} = 0.039;
is_@{co} = 0.05; % import share
cN_c@{co} = 0.53 ; % share of NT to the consumption basket 
cT_c@{co} = 1 - cN_c@{co} ;
% Utilisation parametres
u_ss_@{co} = 1;
elasticity_deltaprime_@{co} = 0.15; % + kk/100 ;
@#endfor
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           DETERMINATION OF THE STEADY-STATE                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Using fsolve for some parameters
x0 = [0.5 ; 0.5 ; 0.5 ; 0.5 ; 0.5 ; 0.5 ; 0.5] ;

%options = optimset('Diagnostics', 'on');
[x,fval, exitflag] = evalin('base', 'fsolve(@myfunCapitalServices3,x0)') ; % Call solver

%if exitflag ~= 1
%    error('NO convergence in fsolve');
%end

% Values solved in fsolve
for i = 1:length(x) ;
       if abs(imag(x(i))) > 1e-8 ;
           disp('Careful: Significant imaginary part') ; 
           chisOnes(i,1) = 0 ;
       end
       chisOnes(i,1) = 1 ;
end
           
if chisOnes == ones(length(x), 1) ;
    disp('GOOD: The imaginary part is almost non-existant') ;
    yT_y1 = real(x(1)) ;
    y1_S1 = real(x(2)) ;
    yT1_ST1 = real(x(3)) ;
    pNyN1_SN1 = real(x(6)) ;
end

for i = 1:length(x)
    if real(x(i)) < 0
        error('Negative Value at ss - from fsolve') ;
    end
end

yN_y1 = 1 - yT_y1 ;

omega_1 = ((yT_y1/is_1 - 1)^(1/theta_1))/(1+(yT_y1/is_1 - 1)^(1/theta_1)) ;
qa1_ss = ((omega_1*(yT_y1 - is_1)^((theta_1-1)/theta_1) + (1-omega_1)*is_1^((theta_1-1)/theta_1))^(theta_1/(theta_1-1)))/yT_y1 ;
qb1_ss = qa1_ss ;

mi1_phi1 = (omega_1*(yT_y1 - is_1)^(-1/theta_1))*((omega_1*(yT_y1-is_1)^((theta_1-1)/theta_1) + (1-omega_1)*is_1^((theta_1-1)/theta_1))^(1/(theta_1-1))) ;

%Investment share
i1_y1 = delta_1/y1_S1 ;
cT_y1 = yT_y1 - i1_y1 ; 

deltaprime_1 = mi1_phi1*alphaT_1*(yT1_ST1^((epsilon_1-1)/epsilon_1))*((yT_y1*y1_S1)^(1/epsilon_1))/qa1_ss;
deltaprime_1_1 = alphaN_1*(pNyN1_SN1^((epsilon_1-1)/epsilon_1))*((yN_y1*y1_S1)^(1/epsilon_1)) ;
deltaprimeDiff = abs(deltaprime_1 - deltaprime_1_1) ;
if abs(deltaprime_1 - deltaprime_1_1) > 1e-5
    error('Check deltaprime')
    disp(['ATTENTION: Check deltaprime, diff is ',  num2str(deltaprimeDiff)]);
end

deltadoubleprime_1 = elasticity_deltaprime_1*deltaprime_1/u_ss_1 ;
nT1_yT1 = ((1/qa1_ss)^(1/(1-alphaT_1)))*((yT1_ST1)^(alphaT_1/(1-alphaT_1))) ;
tau1_mi1 = (1-alphaT_1)/(qa1_ss*nT1_yT1) ;
pNyN1_nN1 = tau1_mi1*mi1_phi1/(1-alphaN_1) ;
p1N_ss = ((pNyN1_SN1)^alphaN_1)*(pNyN1_nN1^(1-alphaN_1)) ;
ksi1_phi1 = p1N_ss ;
nN1_yN1 = (pNyN1_nN1/p1N_ss)^(-1) ;
yN1_SN1 = pNyN1_SN1/p1N_ss ;
nN1_yN1_1 = (yN1_SN1)^(alphaN_1/(1-alphaN_1)) ;
n1_y1  = nT1_yT1*yT_y1 + nN1_yN1*yN_y1/p1N_ss ;
n_ss_1 = 0.33 ;

if abs(nN1_yN1 - nN1_yN1_1) > 1e-7 ;
    error('Check nN1_yN1') 
end;

% Determining omegaT based on ouptut shares. 
omegaT_1 = ((cT_y1/yN_y1)^(1/rho_1))/(p1N_ss^((rho_1-1)/rho_1) + (cT_y1/yN_y1)^(1/rho_1)) ;

% Consumption share
c1_y1 = (omegaT_1*cT_y1^((rho_1-1)/rho_1) + (1-omegaT_1)*(yN_y1/p1N_ss)^((rho_1-1)/rho_1))^(rho_1/(rho_1-1)) ;
% Price of consumption
P1_ss = (cT_y1 + yN_y1)/c1_y1 ;

% Getting PSI right
ksi1_pi1 = ((1-omegaT_1)*(yN_y1/p1N_ss)^(-1/rho_1))*((omegaT_1*cT_y1^((rho_1-1)/rho_1) + (1-omegaT_1)*(yN_y1/p1N_ss)^((rho_1-1)/rho_1))^(1/(rho_1-1))) ;
PSI_1 = (1-alphaN_1)*(1/nN1_yN1)*(1/n_ss_1^(ni_1-1))*(1/ni_1)*ksi1_pi1 ;

P1_ss1 =(1/omegaT_1)*(cT_y1^(1/rho_1))*((omegaT_1*(cT_y1^((rho_1-1)/rho_1)) + (1-omegaT_1)*(yN_y1/p1N_ss)^((rho_1-1)/rho_1))^(-1/(rho_1-1))) ;

disp('         ') % I do that so that it leaves a gap

disp('----------- Steady-state -----------')
% Exact values at ss

% Outputs
y1_ss = n_ss_1*(1/n1_y1) ;
yT1_ss = yT_y1*y1_ss ;
yN1_ss = yN_y1*y1_ss/p1N_ss ;
diseq_out = y1_ss - (yT1_ss + p1N_ss*yN1_ss) ;
if abs(diseq_out) < 1e-08 ;
    disp('Output definition verified: Y = yT + pN*yN')
else
    disp(['ERROR in Output Equations - Output definition NOT verified: Y =/ yT + pN*yN. Difference is ', num2str(diseq_out)])
end

% Labour
nN1_ss = nN1_yN1*yN1_ss ;
nT1_ss = nT1_yT1*yT1_ss ;
diseq_labour = n_ss_1 - (nT1_ss + nN1_ss) ;
if abs(diseq_labour) < 1e-08 ;
    disp('Labour equilibrium verified')
else
    disp(['ERROR in Labour Equations - Labour equilibrium NOT verified. Difference is ', num2str(diseq_labour)]) 
end

% Capital services
S1_ss = (1/y1_S1)*y1_ss ;
ST1_ss = (1/yT1_ST1)*yT1_ss ;
SN1_ss = (1/yN1_SN1)*yN1_ss ;
diseq_KS = S1_ss - (ST1_ss^((epsilon_1-1)/epsilon_1) + SN1_ss^((epsilon_1-1)/epsilon_1))^(epsilon_1/(epsilon_1-1)) ;
if abs(diseq_KS) <1e-06 ;
    disp('Capital Services equilibrium verified')
else
    disp(['ERROR in Capital Services Equations - Capital Services equilibrium NOT verified. Difference is ', num2str(diseq_KS)])
end

% Capital
k1_ss = S1_ss/u_ss_1;

% Consumption and investment
i1_ss = i1_y1*y1_ss ;
c1_ss = c1_y1*y1_ss ;
diseq = y1_ss - (i1_ss + P1_ss*c1_ss) ;
if abs(diseq) < 1e-08 ;
    disp('Equilibrium verified: Y = I + p*C')
else
    disp(['ERROR - Equilibrium NOT verified: P*C + I =/ Y. Difference is ', num2str(diseq)]);
end

% Consumption shares
cT1_ss = cT_y1*y1_ss ;
cN1_ss = yN_y1*y1_ss/p1N_ss ;
%cN1_ss = yN1_ss ;  
diseq_cons = P1_ss*c1_ss - (cT1_ss + p1N_ss*cN1_ss) ;
if abs(diseq_cons) < 1e-08 ;
    disp('Consumption Equilibrium verified: PC = cT + pN*cN')
else
    disp(['ERROR - Consumption Equilibrium NOT verified: PC =/ cT + pN*cN. Difference is ', num2str(diseq_cons)])
end

diseq_NT = yN1_ss - cN1_ss ;
if abs(diseq_NT) < 1e-08 ;
    disp('Equilibrium of NT verified: cN = yN')
else
    disp(['ERROR in NT goods - Equilibrium NOT verified cN =/ yN. Difference is ', num2str(diseq_NT)])
end

% Tradeables
dise_T = yT1_ss - (i1_ss + cT1_ss) ;
if abs(diseq) < 1e-08 ;
    disp('Equilibrium in T verified: yT = cT + I')
else
    disp(['ERROR in T-goods - Equilibrium NOT verified: cT + I =/ yT. Difference is ', num2str(diseq_T)]);
end

yT1_I_ss = yT1_ss/qa1_ss ;
yT2_I_ss = yT1_ss/qb1_ss ;
P_diff = P1_ss - P1_ss1 ;
disp(['The difference in price of final consumption is ', num2str(P_diff)])

if abs(diseq_KS) > 1e-05 | abs(diseq_out)>1e-08 | abs(diseq_labour)>1e-08 | abs(diseq)>1e-08| abs(diseq_cons)>1e-08 | abs(diseq_NT)>1e-08 | abs(P_diff)>1e-08
    error('Check your steady-state')
end; 

% Intermediate goods at ss
a1_ss = (yT_y1 - is_1)*y1_ss/qa1_ss ;
b1_ss = is_1*y1_ss/qa1_ss;

disp('         ') % I do that so that it leaves a gap

% Ratios at Steady State
cT_c1_1 = cT1_ss/(P1_ss*c1_ss);
cN_c1_1 = p1N_ss*cN1_ss/(P1_ss*c1_ss);

if abs(cT_c1_1 - cT_c1) > 1e-06
    error('Something is wrong with consumption shares cT_c1 =/ cT_c1_1')
end

if abs(cN_c1_1 - cN_c1) > 1e-06
    error('Something is wrong with consumption shares: cN_c1 =/ cN_c1_1')
end

if abs(cT_c1_1 - cT_c1) > 1e-08
    disp(['Just a warning about consumption shares: cT_c1 - cT_c1_1 = ', num2str(abs(cT_c1_1-cT_c1))])
end

if abs(cN_c1_1 - cN_c1) > 1e-08
disp(['Just a warning about consumption shares: cN_c1 - cN_c1_1 = ', num2str(abs(cN_c1_1-cN_c1))])
end
disp('   ');

disp('Ratios at steady state')
disp(['Consumption share, PC/Y:', num2str(P1_ss*c1_ss/y1_ss)])
disp(['Investment share, I/Y:', num2str(i1_ss/y1_ss)])
disp(['Share of Tradeables to consumption Basket:', num2str(cT_c1)])
disp(['Share of Non-Tradeables to consumption Basket:', num2str(cN_c1)])
disp(['Output-Capital Services Ratio:', num2str(y1_S1)])

disp('      ')

disp('Calibrated values')
disp(['Elasticity of substitution between home and foreign goods, theta:', num2str(theta_1)])
disp(['Elasticity of substitution between tradeable and non-tradeable goods, rho:', num2str(rho_1)])
disp('Other')
disp(['Share of Non-Tradeables in production:', num2str(yN_y1)])
disp(['Price of the Non-Tradeable good at ss is ', num2str(p1N_ss)])
disp(['Price of Final Consumption at ss is ', num2str(P1_ss)])
h = sprintf('Labour at steady-state is %s with PSI equal to %s.', num2str(n_ss_1), num2str(PSI_1));
disp(h)

% The steady state in a vector

MUc1_ss = (c1_ss - PSI_1*n_ss_1^ni_1)^(gamma_1-1) ; % should I rise to a power? in that case, whether MU is positive or negative depends on gamma. with gamma=-1, this will always be positive. does that make sense?
MUc1_ss1 = (c1_ss - PSI_1*n_ss_1^ni_1) ;
MUn1_ss = ni_1*PSI_1*(n_ss_1^(ni_1-1))*MUc1_ss ;
MUn1_ss1 = ni_1*PSI_1*(n_ss_1^(ni_1-1))*MUc1_ss1 ;
p1N_ss1 = ((1-omegaT_1)/omegaT_1)*(cN1_ss/cT1_ss)^(-1/rho_1) ;
jrT_NT_fixCapacity_ss = [y1_ss yT1_ss yN1_ss n_ss_1 nN1_ss nT1_ss S1_ss ST1_ss SN1_ss i1_ss c1_ss cT1_ss cN1_ss yT1_I_ss qa1_ss P1_ss p1N_ss MUc1_ss MUn1_ss MUc1_ss1 MUn1_ss1]' ;
jrT_NT_fixCapacity_ss_names = ['y1_ss   ' 
                            'yT1_ss  ' 
                            'yN1_ss  ' 
                            'n_ss_1  ' 
                            'nN1_ss  ' 
                            'nT1_ss  '
                            'k1_ss   '
                            'S1_ss   '
                            'ST1_ss  '
                            'SN1_ss  ' 
                            'i1_ss   '
                            'c1_ss   '
                            'cT1_ss  '
                            'cN1_ss  '
                            'yT1_I_ss'
                            'qa1_ss  '
                            'P1_ss   '
                            'p1N_ss  '
                            'uT_ss_1 '
                            'uN_ss_1 '
                            'MUc1_ss '
                            'MUn1_ss '
                            'MUc1_ss1'
                            'MUn1_ss1'] ;

for i = 1:length(jrT_NT_fixCapacity_ss-2) ;
    if jrT_NT_fixCapacity_ss(i) < 0 ;
        error(['ERROR: Negative quantities at ss: ', jrT_NT_fixCapacity_ss_names(i, :) ])
        Vones1(i) = 0;
    end
Vones1(i,1) = 1;
end

disp('       ');

if Vones1 == ones(length(jrT_NT_fixCapacity_ss), 1) ;
      disp('GOOD! All quantities and prices are positive at ss');
end;

% Prints a warning if (c-psi*n^ni) < 0. in that case you will have complex eigenvalues nomizw. 
if MUc1_ss1 < 0 ;
    disp('WARNING: (C - psi*N^ni) is negative. You might have comlex eigenvalues')
end ;


% Langrange Multipliers at ss
pi1_ss = MUc1_ss ;
tau1_ss = ni_1*PSI_1*n_ss_1^(ni_1-1)*pi1_ss ;
phi1_ss = pi1_ss*omegaT_1*(cT1_ss^(-1/rho_1))*(omegaT_1*cT1_ss^((rho_1-1)/rho_1) + (1-omegaT_1)*cN1_ss^((rho_1-1)/rho_1))^(1/(rho_1-1));
%lamda1_ss = phi1_ss ; I want to check this, instead of imposing it
mi1_ss = phi1_ss*(omega_1*a1_ss^(-1/theta_1))*(omega_1*a1_ss^((theta_1-1)/theta_1) + (1-omega_1)*b1_ss^((theta_1-1)/theta_1))^(1/(theta_1-1)) ;
mi1_ss1 = phi1_ss*((1-beta_1*(1-delta_1))*qa1_ss/(beta_1*alphaT_1*(y1_S1*yT_y1)^(1/epsilon_1)))*(1/yT1_ST1)^((epsilon_1-1)/epsilon_1) ;
ksi1_ss = pi1_ss*(1-omegaT_1)*(cN1_ss^(-1/rho_1))*(omegaT_1*cT1_ss^((rho_1-1)/rho_1) + (1-omegaT_1)*cN1_ss^((rho_1-1)/rho_1))^(1/(rho_1-1)) ;
ksi1_ss1 = p1N_ss*phi1_ss ;    
zita1_ss = ksi1_ss*alphaN_1*(SN1_ss^(1/epsilon_1))*yN1_SN1*(ST1_ss^((epsilon_1-1)/epsilon_1) + SN1_ss^((epsilon_1-1)/epsilon_1))^(-1/(epsilon_1-1)) ;
zita1_ss1 = mi1_ss*alphaT_1*(yT1_ST1/qa1_ss)*(ST1_ss^(1/epsilon_1))*(ST1_ss^((epsilon_1-1)/epsilon_1) + SN1_ss^((epsilon_1-1)/epsilon_1))^(-1/(epsilon_1-1)) ; ;
tau1_ss1 = ksi1_ss*(1-alphaN_1)/nN1_yN1 ;
tau1_ss3 = mi1_ss*(1-alphaT_1)*yT1_I_ss/nT1_ss ;
tau1_ss2 = mi1_ss*(1-alphaT_1)/(qa1_ss*nT1_yT1) ;
%tau1_ss2 = mi1_ss*(1-alphaT_1)*yT1_I_ss/nT1_ss ;
lamda1_ss1 = beta_1*zita1_ss*u_ss_1/(1-beta_1*(1-delta_1));
lamda1_ss2 = zita1_ss/deltaprime_1;

jrT_NT_fixCapacity_LMss = [pi1_ss tau1_ss phi1_ss mi1_ss ksi1_ss zita1_ss mi1_ss1 ksi1_ss1  zita1_ss1 tau1_ss1 tau1_ss2 lamda1_ss1 lamda1_ss2]';

jrT_NT_fixCapacity_LMss_names = ['pi1_ss      '
                              'tau1_ss     '
                              'tau1_ss1    '
                              'tau1_ss2    '
                              'tau1_ss3    '
                              'phi1_ss     '
                              'mi1_ss      ' 
                              'ksi1_ss     '
                              'zita1_ss    ' 
                              'mi1_ss1     '
                              'ksi1_ss1    '
                              'zita1_ss1   '
                              'lamda1_ss1  '
                              'lamda1_ss2  '];

for j = 1:length(jrT_NT_fixCapacity_LMss) ;
    if jrT_NT_fixCapacity_LMss(j) < 0 ;
        error(['ERROR: Negative LMs at ss: ', jrT_NT_fixCapacity_LMss_names(j, :) ])
        Vones(j) = 0;
    end
Vones(j,1) = 1; % In order to check if ALL the elements are positive
end

disp('       ');

if Vones == ones(length(jrT_NT_fixCapacity_LMss), 1) ;
      disp('GOOD! All LMs are positive at ss');
else
    erro('You have non-positive LMs at ss');
end;


if abs(mi1_ss - mi1_ss1) > 1e-7 ;
    disp(['ATTENTION: Inequality of LMs at ss. Check mi1_ss. Difference is ', num2str(abs(mi1_ss - mi1_ss1))]);
end

if abs(p1N_ss - p1N_ss1) > 1e-7 ;
    disp(['ATTENTION: Inequality of LMs at ss. Check p1N_ss. Difference is ', num2str(abs(p1N_ss - p1N_ss1))]);
end
if abs(ksi1_ss - ksi1_ss1) > 1e-7  ;
    disp(['ATTENTION: Inequality of LMs at ss. Check ksi1_ss. Difference is ', num2str(abs(ksi1_ss - ksi1_ss1))]);
end
if abs(zita1_ss - zita1_ss1) > 1e-7 ;
    disp(['ATTENTION: Inequality of LMs at ss. Check zita1_ss. Difference is ', num2str(abs(zita1_ss - zita1_ss1))]);
end
if abs(tau1_ss - tau1_ss1) > 1e-7  || abs(tau1_ss - tau1_ss2) > 1e-7 || abs(tau1_ss - tau1_ss3) > 1e-7 ;
    disp(['ATTENTION: Inequality of LMs at ss. Check tau1_ss. Difference is ', num2str(abs(tau1_ss - tau1_ss1))]);
    disp(['ATTENTION: Inequality of LMs at ss. Check tau1_ss. Difference is ', num2str(abs(tau1_ss - tau1_ss2))]);
    disp(['ATTENTION: Inequality of LMs at ss. Check tau1_ss. Difference is ', num2str(abs(tau1_ss - tau1_ss3))]);
end
if abs(lamda1_ss1 - lamda1_ss2) > 1e-7;
    disp(['ATTENTION: Inequality of LMs at ss. Check lamda1_ss. Difference is ', num2str(abs(lamda1_ss1 - lamda1_ss2))]);
end
if abs(lamda1_ss1 - phi1_ss) > 1e-7 || abs(lamda1_ss2 - phi1_ss) > 1e-7;
    disp(['ATTENTION: Inequality of LMs at ss. lamda is not equal to phi. Difference is ', num2str(abs(lamda1_ss1 - phi1_ss))]);
    disp(['ATTENTION: Inequality of LMs at ss. lamda is not equal to phi. Difference is ', num2str(abs(lamda1_ss2 - phi1_ss))]);
end

%if abs(mi1_ss - mi1_ss1) < 1e-4 && abs(ksi1_ss - ksi1_ss1) < 1e-7 && abs(zita1_ss - zita1_ss1) < 1e-4 && abs(tau1_ss - tau1_ss1) < 1e-7  && abs(tau1_ss - tau1_ss2) < 1e-7 && abs(lamda1_ss1 - lamda1_ss2) < 1e-3 && abs(lamda1_ss1 - phi1_ss) < 1e-3 && abs(p1N_ss - p1N_ss1) < 1e-7
%  disp('GOOD! Equality of LMs at ss');
%else 
%  error('Inequality of LMs at ss');  
%end ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Country 2
% Full symmetry
alphaT_2 = alphaT_1 ; 
alphaN_2 = alphaN_1 ;
beta_2 = beta_1 ;
delta_2 = delta_1 ;
gamma_2 = gamma_1 ;
theta_2 = theta_1 ;
epsilon_2 = epsilon_1 ;
rho_2 = rho_1 ;
ni_2 = ni_1 ;
qa2_ss = qa1_ss ;
qb2_ss = qb1_ss ;
yT_y2 = yT_y1 ;
yN_y2 = yN_y1 ;
is_2 = is_1 ;
omega_2 = omega_1 ;
omegaT_2 = omegaT_1 ;
n_ss_2 = n_ss_1 ;
u_ss_2 = u_ss_1 ;
yT2_ST2 = yT1_ST1 ;
y2_S2 = y1_S1 ;
pNyN2_SN2 = pNyN1_SN1;

p2T_ss = qa2_ss*yT_y2*((1-omega_2)*is_2^((theta_2-1)/theta_2) + omega_2*((yT_y2-is_2)^((theta_2-1)/theta_2)))^(-theta_2/(theta_2-1)) ;

mi2_phi2 = (omega_2*(yT_y2 - is_2)^(-1/theta_2))*((omega_2*(yT_y2-is_2)^((theta_2-1)/theta_2) + (1-omega_2)*is_2^((theta_2-1)/theta_2))^(1/(theta_2-1))) ;
i2_y2 = delta_2/y2_S2 ;
cT_y2 = yT_y2 - i2_y2 ;
deltaprime_2 = (mi2_phi2*alphaT_2*yT2_ST2/qa2_ss)*((yT2_ST2/(y2_S2*yT_y2))^(-1/epsilon_1))*((yT_y2*y2_S2/yT2_ST2)^((epsilon_2-1)/epsilon_2) + (yN_y2*y2_S2/pNyN2_SN2)^((epsilon_2-1)/epsilon_2))^(-1/(epsilon_2-1));
deltaprime_2_1 = (alphaN_2*pNyN2_SN2)*((pNyN2_SN2/(y2_S2*yN_y2))^(-1/epsilon_2))*(((yT_y2*y2_S2/yT2_ST2)^((epsilon_2-1)/epsilon_2) + (yN_y2*y2_S2/pNyN2_SN2)^((epsilon_2-1)/epsilon_2))^(-1/epsilon_2));
deltadoubleprime_2 = elasticity_deltaprime_2*deltaprime_2/u_ss_2 ;
% with knowledge of omegaT I can estimate p2N_ss directly
p2N_ss = (((1-omegaT_2)/omegaT_2)^(rho_2/(rho_2-1)))*((cT_y2/yN_y2)^(1/(rho_2-1)))*p2T_ss^(rho_2/(rho_2-1)) ; 
ksi2_phi2 = p2N_ss ;
yN2_SN2 = pNyN2_SN2/p2N_ss;
c2_y2 = (omegaT_2*(cT_y2/p2T_ss)^((rho_2-1)/rho_2) + (1-omegaT_2)*(yN_y2/p2N_ss)^((rho_2-1)/rho_2))^(rho_2/(rho_2-1)) ;
P2_ss = (cT_y2 + yN_y2)/c2_y2 ;
cT_c2 = cT_y2/(P2_ss*c2_y2) ;
cN_c2 = 1-cT_c2 ;
nN2_yN2 = yN2_SN2^(alphaN_2/(1-alphaN_2)) ;
nT2_yT2 = ((1/qa2_ss)^(1/(1-alphaT_2)))*((yT2_ST2)^(alphaT_2/(1-alphaT_2))) ;
n2_y2  = nT2_yT2*yT_y2/p2T_ss + nN2_yN2*yN_y2/p2N_ss ;
% Getting PSI right
ksi2_pi2 = ((1-omegaT_2)*(yN_y2/p2N_ss)^(-1/rho_2))*((omegaT_2*(cT_y2/p2T_ss)^((rho_2-1)/rho_2) + (1-omegaT_2)*(yN_y2/p2N_ss)^((rho_2-1)/rho_2))^(1/(rho_2-1))) ;
PSI_2 = (1-alphaN_2)*ksi2_pi2/(nN2_yN2*ni_2*n_ss_2^(ni_2-1)) ;
P2_ss1 =(1/omegaT_2)*((cT_y2/p2T_ss)^(1/rho_2))*((omegaT_2*((cT_y2/p2T_ss)^((rho_2-1)/rho_2)) + (1-omegaT_2)*(yN_y2/p2N_ss)^((rho_2-1)/rho_2))^(-1/(rho_2-1))) ;

disp('      ');

if abs(mi2_phi2-mi1_phi1)< 1e-10 && abs(yT2_ST2-yT1_ST1)< 1e-10 && abs(p2N_ss-p1N_ss)< 1e-10 &&abs(ksi2_phi2-ksi1_phi1)< 1e-10 && abs(yN2_SN2-yN1_SN1)< 1e-10 && abs(y2_S2-y1_S1)< 1e-10 && abs(i2_y2-i1_y1)< 1e-10 && abs(cT_y2-cT_y1)< 1e-10 && abs(omegaT_2-omegaT_1)< 1e-10 && abs(c2_y2-c1_y1)< 1e-10 && abs(P2_ss-P1_ss)< 1e-10 && abs(nN2_yN2-nN1_yN1)< 1e-10 && abs(nT2_yT2-nT1_yT1)< 1e-10 && abs(n2_y2-n1_y1)< 1e-10 && abs(n_ss_2 - n_ss_1)< 1e-10 && abs(ksi2_pi2-ksi1_pi1)< 1e-10 && abs(PSI_2-PSI_1)< 1e-10 && abs(P1_ss1-P2_ss1)< 1e-10 ;
    disp('Fully symmetric Steady-State')
else
    error('ERROR: NOT-symmetric steady-state')
end;

disp(['Price of the foreign final tradeable good is ', num2str(p2T_ss)]) ;

disp('       ---------------------------      ');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       Parametrisation of the exogenous processes                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Calibration following Corsetti et al
%/*
A1T_11 = 0.82  ;   
A1T_12 = -0.06 ; 
A1T_13 = 0.1   ; 
A1T_14 = 0.24  ;
A2T_11 = A1T_12;
A2T_12 = A1T_11;
A2T_13 = A1T_14;
A2T_14 = A1T_13;
A1N_11 = -0.02 ;
A1N_12 = 0.02  ;
A1N_13 = 0.96  ; 
A1N_14 = 0.01  ;
A2N_11 = A1N_12;
A2N_12 = A1N_11;
A2N_13 = A1N_14;
A2N_14 = A1N_13;
%*/

% The VAR(1) matrix - persistence of the exogenous VAR process
Apersistence = [A1T_11 A1T_12 A1T_13 A1T_14 ; A2T_11 A2T_12 A2T_13 A2T_14 ; A1N_11 A1N_12 A1N_13 A1N_14 ; A2N_11 A2N_12 A2N_13 A2N_14];

% Variance-Covariance Matrix - Corsetti et all
@#for co in countries
sigmaT_@{co} = sqrt(0.047/100); 
sigmaN_@{co}= sqrt(0.009/100);
%seT_@{co} = sqrt(0.047/100) ;
%seN_@{co} = sqrt(0.009/100) ;
@#endfor

% Correlation of shocks - Corsetti et al
correT = 0.022/100 ; %/4.7 ; % tradeables, across countries
correN = -0.001/100; %/0.9 ; % non-tradeables, across countries
correTN1 = 0.009/100; %/(sqrt(4.7)*sqrt(0.9)) ; % T and NT, country one
correTN2 = 0.009/100; %/(sqrt(4.7)*sqrt(0.9)) ; % T and NT, country two
correTNa = 0.004/100;%/(sqrt(4.7)*sqrt(0.9)) ; % T and NT, across countries --> corr(eT1,eN2)
correTNb = 0.004/100;% /(sqrt(4.7)*sqrt(0.9)) ; % T and NT, across countries --> corr(eT2, eN1)

% Correlation of shocks - No positive correlations whatsoever, keep the variances as in Corsetti et al
%correT = 0 ; %/4.7 ; % tradeables, across countries
%correN = 0; %/0.9 ; % non-tradeables, across countries
%correTN1 = 0; %/(sqrt(4.7)*sqrt(0.9)) ; % T and NT, country one
%correTN2 = 0; %/(sqrt(4.7)*sqrt(0.9)) ; % T and NT, country two
%correTNa = 0;%/(sqrt(4.7)*sqrt(0.9)) ; % T and NT, across countries --> corr(eT1,eN2)
%correTNb = 0;% /(sqrt(4.7)*sqrt(0.9)) ; % T and NT, across countries --> corr(eT2, eN1)

% No international Correlation of shocks - Maintain the within country correlations
%correT = 0 ; %/4.7 ; % tradeables, across countries
%correN = 0; %/0.9 ; % non-tradeables, across countries
%correTN1 = 0.009/100; %/(sqrt(4.7)*sqrt(0.9)) ; % T and NT, country one
%correTN2 = 0.009/100; %/(sqrt(4.7)*sqrt(0.9)) ; % T and NT, country two
%correTNa = 0;%/(sqrt(4.7)*sqrt(0.9)) ; % T and NT, across countries --> corr(eT1,eN2)
%correTNb = 0;% /(sqrt(4.7)*sqrt(0.9)) ; % T and NT, across countries --> corr(eT2, eN1)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           THE MODEL NON-LINEARISED                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            

%-------------   First-order condtions  ---------------- 
 
model;
exp(pi1) = (exp(c1) - PSI_1*exp(n1)^ni_1)^(gamma_1-1) ;
exp(pi2) = (exp(c2) - PSI_2*exp(n2)^ni_2)^(gamma_2-1) ;

exp(tau1) = exp(pi1)*ni_1*PSI_1*exp(n1)^(ni_1-1) ;
exp(tau2) = exp(pi2)*ni_2*PSI_2*exp(n2)^(ni_2-1) ;
% new variable: the depreciation function, deltaU1 and deltaU2
% deltaprime the PARAMETER is delta1, the coefficient
exp(deltaU1) = delta_1 + deltaprime_1*(exp(u1) - u_ss_1) + deltadoubleprime_1*(exp(u1)-u_ss_1)^2/2; 
exp(deltaU2) = delta_2 + deltaprime_2*(exp(u2) - u_ss_2) + deltadoubleprime_2*(exp(u2)-u_ss_2)^2/2;
exp(deltaprimeU1) = deltaprime_1 + deltadoubleprime_1*(exp(u1) - u_ss_1);
exp(deltaprimeU2) = deltaprime_2 + deltadoubleprime_2*(exp(u2) - u_ss_2);

exp(lamda1) = beta_1*((1-exp(deltaU1))*exp(lamda1(+1)) + exp(zita1(+1))*exp(u1(+1)));
exp(lamda2) = beta_2*((1-exp(deltaU2))*exp(lamda2(+1)) + exp(zita2(+1))*exp(u2(+1)));

% new variable for investment adjustment cost function: sI. Only linearised, not log-linearised
%exp(sI1) = kappa_1*(exp(i1)/exp(i1(-1)) - 1)^2/2;
%exp(sI2) = kappa_2*(exp(i2)/exp(i2(-1)) - 1)^2/2;
sI1 = kappa_1*(exp(i1)/exp(i1(-1)) - 1)^2/2;
sI2 = kappa_2*(exp(i2)/exp(i2(-1)) - 1)^2/2;

%exp(sprimeI1) = kappa_1*(exp(i1)/exp(i1(-1)) - 1);
%exp(sprimeI2) = kappa_2*(exp(i2)/exp(i2(-1)) - 1);
sprimeI1 = kappa_1*(exp(i1)/exp(i1(-1)) - 1);
sprimeI2 = kappa_2*(exp(i2)/exp(i2(-1)) - 1);

%exp(phi1) = exp(lamda1)*(1 - exp(sI1) - exp(i1)*exp(sprimeI1)/exp(i1(-1))) + beta_1*exp(lamda1(+1))*exp(sprimeI1)*(exp(i1(+1))/exp(i1))^2;
%exp(phi2) = exp(lamda2)*(1 - exp(sI2) - exp(i2)*exp(sprimeI2)/exp(i2(-1))) + beta_2*exp(lamda2(+1))*exp(sprimeI2)*(exp(i2(+1))/exp(i2))^2;
exp(phi1) = exp(lamda1)*(1 - sI1 - exp(i1)*sprimeI1/exp(i1(-1))) + beta_1*exp(lamda1(+1))*sprimeI1*(exp(i1(+1))/exp(i1))^2;
exp(phi2) = exp(lamda2)*(1 - sI2 - exp(i2)*sprimeI2/exp(i2(-1))) + beta_2*exp(lamda2(+1))*sprimeI2*(exp(i2(+1))/exp(i2))^2;

exp(zita1) = exp(lamda1)*exp(deltaprimeU1)*exp(u1) ;
exp(zita2) = exp(lamda2)*exp(deltaprimeU2)*exp(u2) ;

% Consumption shares
exp(phi1) = exp(pi1)*omegaT_1*exp(cT1)^(-1/rho_1)*(omegaT_1*exp(cT1)^((rho_1-1)/rho_1) + (1-omegaT_1)*exp(cN1)^((rho_1-1)/rho_1))^(1/(rho_1-1));
exp(phi2) = exp(pi2)*omegaT_2*exp(cT2)^(-1/rho_2)*(omegaT_2*exp(cT2)^((rho_2-1)/rho_2) + (1-omegaT_2)*exp(cN2)^((rho_2-1)/rho_2))^(1/(rho_2-1));

exp(ksi1) = exp(pi1)*(1-omegaT_1)*exp(cN1)^(-1/rho_1)*(omegaT_1*exp(cT1)^((rho_1-1)/rho_1) + (1-omegaT_1)*exp(cN1)^((rho_1-1)/rho_1))^(1/(rho_1-1));
exp(ksi2) = exp(pi2)*(1-omegaT_2)*exp(cN2)^(-1/rho_2)*(omegaT_2*exp(cT2)^((rho_2-1)/rho_2) + (1-omegaT_2)*exp(cN2)^((rho_2-1)/rho_2))^(1/(rho_2-1));

% Marginal products
alphaT_1*exp(mi1)*exp(yT_I1)/exp(ST1) = exp(zita1)*exp(ST1)^(-1/epsilon_1)*exp(S1)^(1/epsilon_1);
alphaT_2*exp(mi2)*exp(yT_I2)/exp(ST2) = exp(zita2)*exp(ST2)^(-1/epsilon_2)*exp(S2)^(1/epsilon_2);
alphaN_1*exp(ksi1)*exp(yN1)/exp(SN1) = exp(zita1)*exp(SN1)^(-1/epsilon_1)*exp(S1)^(1/epsilon_1);
alphaN_2*exp(ksi2)*exp(yN2)/exp(SN2) = exp(zita2)*exp(SN2)^(-1/epsilon_2)*exp(S2)^(1/epsilon_2);
exp(tau1) = exp(mi1)*(1-alphaT_1)*exp(yT_I1)/exp(nT1);
exp(tau2) = exp(mi2)*(1-alphaT_2)*exp(yT_I2)/exp(nT2);
exp(tau1) = exp(ksi1)*(1-alphaN_1)*exp(yN1)/exp(nN1);
exp(tau2) = exp(ksi2)*(1-alphaN_2)*exp(yN2)/exp(nN2);


% Stochastic processes 
zT1 = A1T_11*zT1(-1) + A1T_12*zT2(-1) + A1T_13*zN1(-1) + A1T_14*zN2(-1) + eT1(-@{q}) ;
zT2 = A2T_11*zT1(-1) + A2T_12*zT2(-1) + A2T_13*zN1(-1) + A2T_14*zN2(-1) + eT2(-@{q}) ;
zN1 = A1N_11*zT1(-1) + A1N_12*zT2(-1) + A1N_13*zN1(-1) + A1N_14*zN2(-1) + eN1(-@{q}) ;
zN2 = A2N_11*zT1(-1) + A2N_12*zT2(-1) + A2N_13*zN1(-1) + A2N_14*zN2(-1) + eN2(-@{q}) ;

%exp(zT1) = A1T_11*exp(zT1(-1)) + A1T_12*exp(zT2(-1)) + A1T_13*exp(zN1(-1)) + A1T_14*exp(zN2(-1)) + eT1(-@{q}) ;
%exp(zT2) = A2T_11*exp(zT1(-1)) + A2T_12*exp(zT2(-1)) + A2T_13*exp(zN1(-1)) + A2T_14*exp(zN2(-1)) + eT2(-@{q}) ;
%exp(zN1) = A1N_11*exp(zT1(-1)) + A1N_12*exp(zT2(-1)) + A1N_13*exp(zN1(-1)) + A1N_14*exp(zN2(-1)) + eN1(-@{q}) ;
%exp(zN2) = A2N_11*exp(zT1(-1)) + A2N_12*exp(zT2(-1)) + A2N_13*exp(zN1(-1)) + A2N_14*exp(zN2(-1)) + eN2(-@{q}) ;

%zT1 = exp(A1T_11)*zT1(-1)*exp(A1T_12)*zT2(-1)*exp(A1T_13)*zN1(-1)*exp(A1T_14)*zN2(-1)*exp(eT1(-@{q})) ;
%zT2 = exp(A2T_11)*zT1(-1)*exp(A2T_12)*zT2(-1)*exp(A2T_13)*zN1(-1)*exp(A2T_14)*zN2(-1)*exp(eT2(-@{q})) ;
%zN1 = exp(A1N_11)*zT1(-1)*exp(A1N_12)*zT2(-1)*exp(A1N_13)*zN1(-1)*exp(A1N_14)*zN2(-1)*exp(eN1(-@{q})) ;
%zN2 = exp(A2N_11)*zT1(-1)*exp(A2N_12)*zT2(-1)*exp(A2N_13)*zN1(-1)*exp(A2N_14)*zN2(-1)*exp(eN2(-@{q})) ;


% FOCs w.r.t intermediate goods. Note that when I do not log-linearise the equations the weights the social planner assigns to countries matter.
% As I assume equal weighting these cancel out but it would not be the case in general
exp(mi1) = exp(phi1)*omega_1*exp(a1)^(-1/theta_1)*exp(yT1)^(1/theta_1);
exp(mi1) = exp(phi2)*(1-omega_2)*exp(a2)^(-1/theta_2)*exp(yT2)^(1/theta_2);
exp(mi2) = exp(phi1)*(1-omega_2)*exp(b1)^(-1/theta_1)*exp(yT1)^(1/theta_1);
exp(mi2) = exp(phi2)*omega_2*exp(b2)^(-1/theta_2)*exp(yT2)^(1/theta_2);

%%%%%%%%%%%%%%%%%%%%%%%% Resource Constraints

% Final goods prodution
exp(yT1) = (omega_1*exp(a1)^((theta_1-1)/theta_1) + (1-omega_1)*exp(b1)^((theta_1-1)/theta_1))^(theta_1/(theta_1-1));
exp(yT2) = (omega_2*exp(b2)^((theta_2-1)/theta_2) + (1-omega_2)*exp(a2)^((theta_2-1)/theta_2))^(theta_2/(theta_2-1));

% Intermediate goods production
exp(yT_I1) = exp(zT1)*exp(ST1)^alphaT_1*exp(nT1)^(1-alphaT_1) ;
exp(yT_I2) = exp(zT2)*exp(ST2)^alphaT_2*exp(nT2)^(1-alphaT_2) ;
exp(yN1)   = exp(zN1)*exp(SN1)^alphaN_1*exp(nN1)^(1-alphaN_1) ;
exp(yN2)   = exp(zN2)*exp(SN2)^alphaN_2*exp(nN2)^(1-alphaN_2) ;

%... with:
exp(yT_I1) = exp(a1) + exp(a2);
exp(yT_I2) = exp(b1) + exp(b2);

% This is aggregate output, expressed in each country's final tradeable good.
exp(y1) = exp(yT1) + exp(pN1)*exp(yN1) ;
exp(y2) = exp(yT2) + (exp(pN2)/exp(p2T))*exp(yN2) ;
exp(y1_real) = exp(yT1) + exp(yN1) ;
exp(y2_real) = exp(yT2) + exp(yN2) ;

%... with:
exp(yN1) = exp(cN1) ;
exp(yN2) = exp(cN2) ;
exp(yT1) = exp(cT1) + exp(i1) ;
exp(yT2) = exp(cT2) + exp(i2) ;

% Aggregate consumption
exp(c1) = (omegaT_1*exp(cT1)^((rho_1-1)/rho_1) + (1-omegaT_1)*exp(cN1)^((rho_1-1)/rho_1))^(rho_1/(rho_1-1));
exp(c2) = (omegaT_2*exp(cT2)^((rho_2-1)/rho_2) + (1-omegaT_2)*exp(cN2)^((rho_2-1)/rho_2))^(rho_2/(rho_2-1));

% Investment
%exp(k1) = (1-exp(deltaU1))*exp(k1(-1)) + exp(i1)*(1-exp(sI1));
%exp(k2) = (1-exp(deltaU2))*exp(k2(-1)) + exp(i2)*(1-exp(sI2));
exp(k1) = (1-exp(deltaU1))*exp(k1(-1)) + exp(i1)*(1-sI1);
exp(k2) = (1-exp(deltaU2))*exp(k2(-1)) + exp(i2)*(1-sI2);


% Capital Services
exp(S1) = (exp(ST1)^((epsilon_1-1)/epsilon_1) + exp(SN1)^((epsilon_1-1)/epsilon_1))^(epsilon_1/(epsilon_1-1)) ;
exp(S2) = (exp(ST2)^((epsilon_2-1)/epsilon_2) + exp(SN2)^((epsilon_2-1)/epsilon_2))^(epsilon_2/(epsilon_2-1)) ;

% Labour
exp(n1) = exp(nT1) + exp(nN1) ;
exp(n2) = exp(nT2) + exp(nN2) ;

% Utilisation - Capital services
exp(S1) = exp(u1)*exp(k1(-1)) ;
exp(S2) = exp(u2)*exp(k2(-1)) ;

% Equilibrium

exp(yT1) = exp(qa1)*exp(a1) + exp(qb1)*exp(b1) ;
exp(p2T)*exp(yT2) = exp(qb2)*exp(b2) + exp(qa2)*exp(a2) ;

exp(P1)*exp(c1) = exp(cT1) + exp(pN1)*exp(cN1) ;
exp(P2)*exp(c2) = exp(p2T)*exp(cT2) + exp(pN2)*exp(cN2) ;

% Prices
exp(p) = ((1-omega_1)/omega_1)*(exp(b1)/exp(a1))^(-1/theta_1) ;
exp(p)= exp(qb1)/exp(qa1) ;
exp(qa1) = exp(qa2) ;
exp(qb1) = exp(qb2) ; 
exp(pN1) = ((1-omegaT_1)/omegaT_1)*(exp(cN1)/exp(cT1))^(-1/rho_1);
exp(pN2)/exp(p2T) = ((1-omegaT_2)/omegaT_2)*(exp(cN2)/exp(cT2))^(-1/rho_2);
% The real exchange rate: the ratio of the price of foreign over home aggregate consumption
exp(rerCONS) = exp(P2)/exp(P1); 
% the ratio of prices of the NT good
exp(rerNT) = exp(pN2)/exp(pN1) ;
% the RER = the ratio of LMs <=> ratio of marginal utilities of consumption
exp(rerLM) = exp(pi2)/exp(pi1) ;

% Other - definitions

% relative consumption
exp(c_rel) = exp(c1)/exp(c2) ;
exp(cN1_cT1) = exp(cN1)/exp(cT1) ;

% net exports: linearised, not log-linearised
%nx1_y1 = is_1*(qa2_ss*(exp(qa2) - 1)) + qa2_ss*(is_1*(exp(a2) - 1)) 
%       - is_1*(qb1_ss*(exp(qb1) - 1)) - qb1_ss*(is_1*(exp(b1) - 1)) ;
%nx1_y1_real1 = is_1*(exp(a2)-1) - is_1*(exp(b1)-1);
%nx1_y1 = (qa1*a2 - qb1*b1)/y1_real;

% relative output
exp(y_rel) = exp(y1_real)/exp(y2_real);
% relative labour
exp(n_rel) = exp(n1)/exp(n2) ;

exp(mi_phi) = exp(mi1)/exp(phi1) ;
end;

initval;
p = log(1);
rerCONS = log(1);
p2T = log(1);
c_rel = log(1);
%nx1_y1 = log(1);
y_rel = log(1);
c1= log(c1_ss);
n1= log(n_ss_1);
k1= log(k1_ss);
y1= log(y1_ss);
i1= log(i1_ss);
zT1= log(1);
zN1= log(1);
qa1= log(qa1_ss);
qb1= log(qa1_ss);
a1= log(a1_ss);
b1= log(b1_ss);
mi1= log(mi1_ss);
lamda1= log(lamda1_ss1);
phi1= log(phi1_ss);
cT1= log(cT1_ss);
cN1= log(cN1_ss);
yT1= log(yT1_ss);
yT_I1= log(yT1_I_ss);
yN1= log(yN1_ss);
pN1= log(p1N_ss);
S1= log(S1_ss);
ST1= log(ST1_ss);
SN1= log(SN1_ss);
nT1= log(nT1_ss);
nN1= log(nN1_ss);
P1= log(P1_ss);
u1= log(u_ss_1);
y1_real= log(yT1_ss + yN1_ss);
ksi1= log(ksi1_ss);
pi1= log(pi1_ss);
zita1= log(zita1_ss);
tau1= log(tau1_ss);
deltaU1= log(delta_1);
deltaprimeU1 = log(deltaprime_1);
sI1= 0; //exp(1); // 0; //1; //log(1);
sprimeI1= 0; //exp(1); // 0; //1; // log(1);

% country 2
c2= log(c1_ss);
n2= log(n_ss_1);
k2= log(k1_ss);
y2= log(y1_ss);
i2= log(i1_ss);
zT2= log(1);
zN2= log(1);
qa2= log(qa1_ss);
qb2= log(qb1_ss);
a2= log(b1_ss);
b2= log(a1_ss);
mi2= log(mi1_ss);
lamda2= log(lamda1_ss1 );
phi2 = log(phi1_ss);
cT2= log(cT1_ss);
cN2= log(cN1_ss);
yT2= log(yT1_ss);
yT_I2= log(yT2_I_ss);
yN2= log(yN1_ss);
pN2= log(p2N_ss );
S2= log(S1_ss); 
ST2= log(ST1_ss);
SN2= log(SN1_ss);
nT2= log(nT1_ss);
nN2= log(nN1_ss);
P2= log(P2_ss);
u2= log(u_ss_2);
y2_real= log(yT1_ss + yN1_ss);
ksi2= log(ksi1_ss);
pi2= log(pi1_ss);
zita2= log(zita1_ss);
tau2= log(tau1_ss);
deltaU2= log(delta_2);
deltaprimeU2= log(deltaprime_2);
sI2= 0; //1; //exp(1); //0; //1; //log(1);
sprimeI2= 0; //1; //exp(1); // 0; //1; //log(1);

rerNT = log(p2N_ss/p1N_ss);
rerLM = log(1);
cN1_cT1 = log(cN1_ss/cT1_ss);
n_rel = log(n_ss_1/n_ss_1);
mi_phi= log(mi1_ss/phi1_ss );

end;

resid(1);
steady;
check ;

shocks ;
var eT1 = sigmaT_1^2 ;
var eT2 = sigmaT_2^2 ;
var eN1 = sigmaN_1^2 ;
var eN2 = sigmaN_2^2 ;
var eT1, eT2 = correT ;
var eT1, eN2 = correTNa ;
var eT1, eN1 = correTN1 ;
var eT2, eN1 = correTNb ;
var eT2, eN2 = correTN2 ;
var eN1, eN2 = correN ;
end ;

stoch_simul(order=1, hp_filter=100, nograph, ar=0, irf=10); % p rerCONS c1 c2 c_rel cN1_cT1 nx1_y1 nx1_y1_real1 y1_real y2_real y1 y2 i1 i2 n1 n2 P1 P2 pN1 pN2 rerNT S1 S2 ST1 SN1 yT_I1 y_rel u1 u2 yT1 yN1 cT1 cN1 nT1 nN1 zT1 zT2 zN1 zN2 n_rel qa1 mi_phi k1 k2;