% BS paper with JR preferenes
% model non-linearised, letting dynare do it
% log-linearisation

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

q= @{q}

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

var p rerCONS p2T c_rel nx1_y1 nx1_y1_real1 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} y@{co}_nipa  X@{co} jrLM@{co} jrLMops@{co}; 

varexo eT@{co} eN@{co} ;
@#endfor
var rerNT rerLM cN1_cT1 n_rel mi_phi nx1_y1_nipa y_rel_nipa mpkT_1 mpkN_1 yN1_nipa yT2_nipa yN2_nipa; % yN1_yT1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% 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} gammaJR_@{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
%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;
gammaJR_@{co} = 0.001;
theta_@{co} = 2; 
epsilon_@{co} = 1.67; //1.635; //1.856; //2.05; //2000000000000;
rho_@{co} = 0.74;
ni_@{co} = 1.4 ;
kappa_@{co} = 0.05; //0; % Difference with before, I set kappa a bit higher to be a bit more realistic maybe
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 ; // 200000; % when I use a much higher value to make it fixed-capacity, it gives me a warning about matrix close to singlular. so I leave it a bit smaller, still enough. 
@#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))) ;

% solving for PSI explicitly (no fsolve)
PSI_1 = (1/n_ss_1^ni_1)*(1-alphaN_1)/nN1_yN1*ksi1_pi1*(ni_1*c1_y1/n1_y1 - (1-alphaN_1)/nN1_yN1*ksi1_pi1*gammaJR_1/(beta_1*(1-gammaJR_1)-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;

%NIPA output
y1_nipa_ss = yT1_ss + p1N_ss*yN1_ss ;
% Sectoral output at constant prices
% yT1 = yT1_nipa because it is the numberaire
yN1_nipa_ss = p1N_ss*yN1_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

jrLM_ss1 = PSI_1*n_ss_1^ni_1*(c1_ss*(1 - PSI_1*n_ss_1^ni_1))^(gamma_1-1)/(-1 + beta_1*(1-gammaJR_1)) ;
MUc1_ss = (c1_ss*(1- PSI_1*n_ss_1^ni_1))^(gamma_1-1) + gammaJR_1*jrLM_ss1 ;
MUn1_ss = ni_1*PSI_1*n_ss_1^(ni_1-1)*c1_ss*(c1_ss - PSI_1*n_ss_1^ni_1*c1_ss)^(gamma_1 - 1) ;
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]' ;
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  '] ;

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;

% Langrange Multipliers at ss
pi1_ss = MUc1_ss ;
tau1_ss = MUn1_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 = PSI_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;
*/
% My estimation of TFP processes - Corsetti et al
A1T_11 = 0.9616  ;   
A1T_12 = -0.0402 ; 
A1T_13 = -0.2303 ; 
A1T_14 =  0.1748 ;
A2T_11 = A1T_12;
A2T_12 = A1T_11;
A2T_13 = A1T_14;
A2T_14 = A1T_13;
A1N_11 = 0.0396;
A1N_12 = -0.0083;
A1N_13 = 0.5595 ; 
A1N_14 = 0.0730 ;
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);
@#endfor
*/

% Variance-Covariance Matrix - My estimations
@#for co in countries
%sigmaT_@{co} = 0.6886*1e-4 ; %BT
%sigmaN_@{co} = 0.2145*1e-4 ; %BT
sigmaT_@{co} = 0.5866*1e-4 ; %corsetti
sigmaN_@{co} = 0.1785*1e-4 ; %corsetti
@#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 - ME, Corsetti et al
correT   = 0.0504*1e-4 ;  % tradeables, across countries
correN   = 0.0363*1e-4 ; % non-tradeables, across countries
correTN1 = 0.1477*1e-4 ; % T and NT, country one
correTN2 = 0.1477*1e-4 ; % T and NT, country two
correTNa = 0.0085*1e-4 ; % T and NT, across countries --> corr(eT1,eN2)
correTNb = 0.0085*1e-4 ; % T and NT, across countries --> corr(eT2, eN1)


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

%-------------   First-order condtions  ---------------- 
 
model;
exp(jrLMops1) = -jrLM1;
exp(jrLMops2) = -jrLM2;

exp(pi1) = (exp(c1) - PSI_1*exp(n1)^ni_1*exp(X1))^(gamma_1-1) - gammaJR_1*exp(jrLMops1)*exp(c1)^(gammaJR_1-1)*exp(X1(-1))^(1-gammaJR_1);
exp(pi2) = (exp(c2) - PSI_2*exp(n2)^ni_2*exp(X2))^(gamma_2-1) - gammaJR_2*exp(jrLMops2)*exp(c2)^(gammaJR_2-1)*exp(X2(-1))^(1-gammaJR_2);

% NEW: FOC wrt X
PSI_1*exp(n1)^ni_1*(exp(c1) - PSI_1*exp(n1)^ni_1*exp(X1))^(gamma_1-1) - exp(jrLMops1) = - beta_1*(1-gammaJR_1)*exp(jrLMops1(+1))*exp(c1(+1))^gammaJR_1*exp(X1)^(-gammaJR_1);
PSI_2*exp(n2)^ni_2*(exp(c2) - PSI_2*exp(n2)^ni_2*exp(X2))^(gamma_2-1) - exp(jrLMops2) = - beta_2*(1-gammaJR_2)*exp(jrLMops2(+1))*exp(c2(+1))^gammaJR_2*exp(X2)^(-gammaJR_2);

exp(tau1) = PSI_1*ni_1*exp(n1)^(ni_1-1)*exp(X1)*(exp(c1) - PSI_1*exp(n1)^ni_1*exp(X1))^(gamma_1-1);
exp(tau2) = PSI_2*ni_2*exp(n2)^(ni_2-1)*exp(X2)*(exp(c2) - PSI_2*exp(n2)^ni_2*exp(X2))^(gamma_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(+1)))*exp(lamda1(+1)) + exp(zita1(+1))*exp(u1(+1)));
exp(lamda2) = beta_2*((1-exp(deltaU2(+1)))*exp(lamda2(+1)) + exp(zita2(+1))*exp(u2(+1)));

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

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

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

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

% 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);

% Create variables for the marginal products
mpkT_1 = exp(zita1)/exp(mi1) ;
mpkN_1 = exp(zita1)/exp(ksi1) ;

% 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}) ;

% 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) ;
% Aggregate REAL output
exp(y1_real) = exp(yT1) + exp(yN1) ;
exp(y2_real) = exp(yT2) + exp(yN2) ;
% Aggregate output at CONSTANT prices
exp(y1_nipa) = exp(yT1) + p1N_ss*exp(yN1) ;
exp(y2_nipa) = p2T_ss*exp(yT2) + p2N_ss*exp(yN2) ;
% Sectoral output at constant prices
% yT1 = yT1_nipa because it is the numberaire
exp(yN1_nipa) = p1N_ss*exp(yN1) ;
exp(yT2_nipa) = p2T_ss*exp(yT2) ;
exp(yN2_nipa) = p2N_ss*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-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) ;

exp(X1) = exp(c1)^gammaJR_1*exp(X1(-1))^(1-gammaJR_1);
exp(X2) = exp(c2)^gammaJR_2*exp(X2(-1))^(1-gammaJR_2);

% 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 = (exp(qa1)*exp(a2) - exp(qb1)*exp(b1))/exp(y1_real);
nx1_y1_real1 = (exp(a2) - exp(b1))/exp(y1_real);
nx1_y1_nipa = (qa1_ss*exp(a2) - qb1_ss*exp(b1))/exp(y1_nipa);

% relative output
exp(y_rel) = exp(y1_real)/exp(y2_real);
exp(y_rel_nipa) = exp(y1_nipa)/exp(y2_nipa);
//exp(yN1_yT1) = exp(yN1)/exp(yT1) ;

% relative labour
exp(n_rel) = exp(n1)/exp(n2) ;

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

initval;
jrLM1 = jrLM_ss1;
jrLMops1 = log(-jrLM_ss1);
p = log(1);
rerCONS = log(1);
p2T = log(1);
c_rel = log(1);
y1_nipa = log(y1_nipa_ss);
y_rel = log(1);
c1= log(c1_ss);
X1 = 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_ss1);
deltaU1= log(delta_1);
deltaprimeU1 = log(deltaprime_1);
sI1= 0; 
sprimeI1= 0;

% country 2
jrLM2 = jrLM_ss1;
jrLMops2 = log(-jrLM_ss1);
c2= log(c1_ss);
X2= log(c1_ss);
n2= log(n_ss_1);
k2= log(k1_ss);
y2_nipa = log(y1_nipa_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; 
sprimeI2= 0;

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 );
//yN1_yT1 = log(1);

yN1_nipa = log(p1N_ss*yN1_ss) ;
yT2_nipa = log(p2T_ss*yT1_ss) ;
yN2_nipa = log(p2N_ss*yN1_ss) ;

end;

resid(1);
steady(solve_algo = 0);
check ;

shocks ;
var eT1 = sigmaT_1 ;
var eT2 = sigmaT_2 ;
var eN1 = sigmaN_1 ;
var eN2 = sigmaN_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 y_rel_nipa nx1_y1_nipa c_rel mpkT_1 mpkN_1 ST1 SN1 yT1 yN1_nipa yN1 yT2_nipa yN2_nipa yT2 nT1 nN1 ;
% 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 nx1_y1_nipa y_rel_nipa y1_nipa y2_nipa ;
 
%p rerCONS y_rel_nipa nx1_y1_nipa c_rel; 

disp('     ');

disp(['Note 1: Theta is ', num2str(theta_1)])
disp(['Note 2: Epsilon is ', num2str(epsilon_1)])
disp(['Note 3: Interim period is ', num2str(q)])
disp(['Note 7: Investment adjustment cost parametre is  ', num2str(kappa_1)]);
disp(['Note 4: ni is ', num2str(ni_1)])
disp(['Note 5: Import share is ', num2str(is_1)])
if abs(A1T_11 - 0.82)< 1e-12 && abs(sigmaT_1- sqrt(0.047/100))< 1e-12 && abs(sigmaN_1 - sqrt(0.009/100))< 1e-12 && abs(correT - 0.022/100)< 1e-12 ;
    disp('Note 5: Shock parametrisation as in Corsetti et al')
elseif abs(A1T_11-0.84)< 1e-12 && abs(sigmaT_1 - sqrt(3.76/1000))< 1e-12 && abs(sigmaN_1 - sqrt(0.51/1000))< 1e-12 && abs(correT - 1.59/1000)< 1e-12 ;
    disp('Note 5: Shock parametrisation as in BT')
end;
disp(['Note 6: Elasticity of deltaprime is  ', num2str(deltadoubleprime_1*u_ss_1/deltaprime_1)])
disp(['Note 8: GammaJR is  ', num2str(gammaJR_1)]);

if alphaT_1 == 0.33
    disp(['Note 7: Labour share in production alphaT following BT: ', num2str(alphaT_1)]) ;
end;

/*
% Keeping the IRFs
% I call it "Benchmark 4" as until now I was calling it 3
%eval(['irfsBenchmark4_', num2str(@{q}), ' = oo_.irfs;']); 
%s = ['save irfsBenchmark4_' num2str(@{q}) '.mat, irfsBenchmark4_' num2str(@{q});] ;
%eval(s) ;
%close all

% Keeping the IRFs, epsilon = 1.45
%eval(['irfsBenchmark4Epsilon145_', num2str(@{q}), ' = oo_.irfs;']); 
%s = ['save irfsBenchmark4Epsilon145_' num2str(@{q}) '.mat, irfsBenchmark4Epsilon145_' num2str(@{q});] ;
%eval(s) ;
%close all

% Keeping the IRFs, epsilon = 3
% eval(['irfsBenchmark4Epsilon3_', num2str(@{q}), ' = oo_.irfs;']); 
% s = ['save irfsBenchmark4Epsilon3_' num2str(@{q}) '.mat, irfsBenchmark4Epsilon3_' num2str(@{q});] ;
% eval(s) ;
% close all

% Keeping the IRFs - "Homogeneous capital": I set epsilon > 2m
% I call it "Benchmark 4" as until now I was calling it 3
% eval(['irfsHomogeneousCapital4_', num2str(@{q}), ' = oo_.irfs;']); 
% s = ['save irfsHomogeneousCapital4_' num2str(@{q}) '.mat, irfsHomogeneousCapital4_' num2str(@{q});] ;
% eval(s) ;
% close all

% Keeping the IRFs - "Fixed capacity": I set elasticity_deltaprime = 200000
% I call it "Benchmark 4" as until now I was calling it 3
% eval(['irfsFixedCapacity4_', num2str(@{q}), ' = oo_.irfs;']); 
% s = ['save irfsFixedCapacity4_' num2str(@{q}) '.mat, irfsFixedCapacity4_' num2str(@{q});] ;
% eval(s) ;
% close all

% Keeping the IRFs - "Fixed capacity": I set elasticity_deltaprime > 2m
% I call it "Benchmark 4" as until now I was calling it 3
%eval(['irfsFixedCapacityHomoCapital4_', num2str(@{q}), ' = oo_.irfs;']); 
%s = ['save irfsFixedCapacityHomoCapital4_' num2str(@{q}) '.mat, irfsFixedCapacityHomoCapital4_' num2str(@{q});] ;
%eval(s) ;
%close all
*/

/*
disp(['Note 1: Theta is ', num2str(theta_1)])
disp(['Note 2: Epsilon is ', num2str(epsilon_1)])
disp(['Note 3: Interim period is ', num2str(q)])
disp(['Note 7: Investment adjustment cost parametre is  ', num2str(kappa_1)]);
disp(['Note 4: ni is ', num2str(ni_1)])
disp(['Note 5: Import share is ', num2str(is_1)])
if abs(A1T_11 - 0.82)< 1e-12 && abs(sigmaT_1- sqrt(0.047/100))< 1e-12 && abs(sigmaN_1 - sqrt(0.009/100))< 1e-12 && abs(correT - 0.022/100)< 1e-12 ;
    disp('Note 5: Shock parametrisation as in Corsetti et al')
elseif abs(A1T_11-0.84)< 1e-12 && abs(sigmaT_1 - sqrt(3.76/1000))< 1e-12 && abs(sigmaN_1 - sqrt(0.51/1000))< 1e-12 && abs(correT - 1.59/1000)< 1e-12 ;
    disp('Note 5: Shock parametrisation as in BT')
end;
disp(['Note 6: Elasticity of deltaprime is  ', num2str(deltadoubleprime_1*u_ss_1/deltaprime_1)])
disp(['Note 7: Investment adjustment cost parametre is  ', num2str(kappa_1)]);

if alphaT_1 == 0.33
    disp(['Note 7: Labour share in production alphaT following BT: ', num2str(alphaT_1)]) ;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%% Extra things - Used to provide me with some tables %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Coefficients of the log linearised RER
Arer = - ((gamma_1-1)*(c1_y1)/(c1_y1-PSI_1*n1_y1*n_ss_1^(ni_1-1)))*sqrt(oo_.var(5,5))/sqrt(oo_.var(2,2));
Brer = (ni_1*PSI_1*(gamma_1-1)*n1_y1*n_ss_1^(ni_1-1)/(c1_y1-PSI_1*n1_y1*n_ss_1^(ni_1-1)))*oo_.var(40,5)/(sqrt(oo_.var(2,2))*sqrt(oo_.var(5,5)));
%/*
%Remember that the matrix oo_.var includes the var-cov matrix of the variables you asked statistics on. 
%Thus, the position of these variables need not be the same as of the M_.endo_names
                                         % position of output: 9 
Names = ['Consumption                 '; % position in var-cov matrix:  3
         'Investment                  '; % position in var-cov matrix: 13
         'Hours                       '; % position in var-cov matrix: 15
         'Real Exchange Rate          '; % position in var-cov matrix:  2
         'Terms of Trade              '; % position in var-cov matrix:  1
         'Relative price of NT        '; % position in var-cov matrix:  21
         'NX over GDP (real)          '; % position in var-cov matrix:  8
         'NX over GDP (current prices)'];% position in var-cov matrix:  7

RelSDs = [sqrt(oo_.var(3,3)/oo_.var(9,9)) ; sqrt(oo_.var(13,13)/oo_.var(9,9)) ; sqrt(oo_.var(15,15)/oo_.var(9,9)) ; sqrt(oo_.var(2,2)/oo_.var(9,9)) ; sqrt(oo_.var(1,1)/oo_.var(9,9)) ; sqrt(oo_.var(21,21)/oo_.var(9,9)) ; sqrt(oo_.var(8,8)/oo_.var(9,9)) ; sqrt(oo_.var(7,7)/oo_.var(9,9))];

% Within country co-movements
NBCs = [oo_.var(3,9)/(sqrt(oo_.var(3,3))*sqrt(oo_.var(9,9))) ; oo_.var(15,9)/(sqrt(oo_.var(15,15))*sqrt(oo_.var(9,9))) ; oo_.var(13,9)/(sqrt(oo_.var(13,13))*sqrt(oo_.var(9,9))) ; oo_.var(8,9)/(sqrt(oo_.var(8,8))*sqrt(oo_.var(9,9))) ];
NamesNBCs = ['Consumption                 '; 
             'Hours                       '; 
             'Investment                  '; 
             'NX over GDP (real)          '];

IBCs = [oo_.var(9,10)/(sqrt(oo_.var(9,9))*sqrt(oo_.var(10,10))) ; oo_.var(3,4)/(sqrt(oo_.var(3,3))*sqrt(oo_.var(4,4))) ; oo_.var(15,16)/(sqrt(oo_.var(15,15))*sqrt(oo_.var(16,16))) ; oo_.var(13,14)/(sqrt(oo_.var(14,14))*sqrt(oo_.var(13,13))) ];
NamesIBCs = ['Outputs (real)              ';  
             'Consumptions                '; 
             'Hours                       '; 
             'Investment                  '];

RERcorrs = [oo_.var(2,5)/(sqrt(oo_.var(2,2))*sqrt(oo_.var(5,5))) ; oo_.var(2,1)/(sqrt(oo_.var(2,2))*sqrt(oo_.var(1,1))) ; oo_.var(2,27)/(sqrt(oo_.var(2,2))*sqrt(oo_.var(27,27))); oo_.var(2,8)/(sqrt(oo_.var(2,2))*sqrt(oo_.var(8,8))) ]; 
NamesRERcorrs = ['Rel. Consumption    ';  
                 'Terms of Trade      '; 
                 'Rel. Output         '; 
                 'NX over GDP (real)  '];

% Value of Loss Function
Loss4 = (RERcorrs(1) - (-0.42))^2 + (RERcorrs(3) - (-0.19))^2 ...
        + (RERcorrs(2) - (0.52))^2 + (RERcorrs(4) - (0.6))^2 ;

Loss3 = (RERcorrs(3) - (-0.19))^2 ...
        + (RERcorrs(2) - (0.52))^2 + (RERcorrs(4) - (0.6))^2 ;

Loss2 = (RERcorrs(2) - (0.52))^2 + (RERcorrs(4) - (0.6))^2 ;

% Loss function: ALL moments that I present in my table

LossAll = (RelSDs(1)-(0.94))^2 + (RelSDs(2)-(4.33))^2 + (RelSDs(3)-(1.19))^2 ...
+ (RelSDs(4)-(3.9))^2 + (RelSDs(5)-(1.68))^2 + (RelSDs(6)-(0.86))^2 + ... 
(NBCs(4)-(-0.48))^2 + (IBCs(1)-(0.68))^2 + (IBCs(2)-(0.60))^2 + (IBCs(3)-(0.54))^2 ...
+ (IBCs(4)-(0.25))^2 + Loss4;   

% Loss function: ALL moments EXCLUDING B-S

LossAllexcludingBS = (RelSDs(1)-(0.94))^2 + (RelSDs(2)-(4.33))^2 + (RelSDs(3)-(1.19))^2 ...
+ (RelSDs(4)-(3.9))^2 + (RelSDs(5)-(1.68))^2 + (RelSDs(6)-(0.86))^2 + ... 
(NBCs(4)-(-0.48))^2 + (IBCs(1)-(0.68))^2 + (IBCs(2)-(0.60))^2 + (IBCs(3)-(0.54))^2 ...
+ (IBCs(4)-(0.25))^2 + Loss3;   

dyntable(['------ Relative Standard Deviations, news are ', num2str(q), ' periods ------'], strvcat('Variable','  Rel. SD'), Names, RelSDs, 10,12,3);
disp('    ');
dyntable(['------ Correlations between GDP and ..., news are ', num2str(q), ' periods ------'], strvcat('Variable','        Correlation'), NamesNBCs, NBCs, 10,12,3);
disp('    ');
dyntable(['------ International Correlations, news are ', num2str(q), ' periods ------'], strvcat('Variable','             Correlation'), NamesIBCs, IBCs, 10,12,3);
disp('    ');
dyntable(['------ Correlation of RER with..., news are ', num2str(q), ' periods ------'], strvcat('Variable','             Correlation'), NamesRERcorrs, RERcorrs, 10,12,3);
disp('    ');

disp('------ The VALUES of the Loss functions are --------')
hl2 = sprintf('Taking into account all moments %s ', num2str(Loss4));
disp(hl2)
hl2 = sprintf('Taking into account 3 moments %s ', num2str(Loss3));
disp(hl2)
hl2 = sprintf('Taking into account 2 moments %s ', num2str(Loss2));
disp(hl2)
hl2 = sprintf('Taking into account EVERYTHING %s ', num2str(LossAll));
disp(hl2)
hl2 = sprintf('Taking into account EVERYTHING but the BS corr %s ', num2str(LossAllexcludingBS));
disp(hl2)

% Keeping the IRFs
%eval(['irfsBenchmark3_', num2str(@{q}), ' = oo_.irfs;']); 
%s = ['save irfsBenchmark3_' num2str(@{q}) '.mat, irfsBenchmark3_' num2str(@{q});] ;
%eval(s) ;
%close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Statistics with NIPA output and NX_Y    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

RelSDs_NIPA = [sqrt(oo_.var(3,3)/oo_.var(47,47)) ; sqrt(oo_.var(13,13)/oo_.var(47,47)) ; sqrt(oo_.var(15,15)/oo_.var(47,47)) ; sqrt(oo_.var(2,2)/oo_.var(47,47)) ; sqrt(oo_.var(1,1)/oo_.var(47,47)) ; sqrt(oo_.var(21,21)/oo_.var(47,47)) ; sqrt(oo_.var(8,8)/oo_.var(47,47)) ; sqrt(oo_.var(7,7)/oo_.var(47,47))];
NBCs_NIPA = [oo_.var(3,46)/(sqrt(oo_.var(3,3))*sqrt(oo_.var(47,47))) ; oo_.var(15,46)/(sqrt(oo_.var(15,15))*sqrt(oo_.var(47,47))) ; oo_.var(13,46)/(sqrt(oo_.var(13,13))*sqrt(oo_.var(47,47))) ; oo_.var(45,47)/(sqrt(oo_.var(45,45))*sqrt(oo_.var(47,47))) ];
IBCs_NIPA = [oo_.var(47,48)/(sqrt(oo_.var(47,47))*sqrt(oo_.var(48,48))) ; oo_.var(3,4)/(sqrt(oo_.var(3,3))*sqrt(oo_.var(4,4))) ; oo_.var(15,16)/(sqrt(oo_.var(15,15))*sqrt(oo_.var(16,16))) ; oo_.var(13,14)/(sqrt(oo_.var(14,14))*sqrt(oo_.var(13,13))) ];
RERcorrs_NIPA = [oo_.var(2,5)/(sqrt(oo_.var(2,2))*sqrt(oo_.var(5,5))) ; oo_.var(2,1)/(sqrt(oo_.var(2,2))*sqrt(oo_.var(1,1))) ; oo_.var(2,46)/(sqrt(oo_.var(2,2))*sqrt(oo_.var(46,46))); oo_.var(2,45)/(sqrt(oo_.var(2,2))*sqrt(oo_.var(45,45))) ]; 

disp('     ');
disp('  --------------------------  NIPA  --------------------------------  ');
dyntable(['------ Relative Standard Deviations, news are ', num2str(q), ' periods ------'], strvcat('Variable','  Rel. SD'), Names, RelSDs_NIPA, 10,12,3);
disp('    ');
dyntable(['------ Correlations between GDP and ..., news are ', num2str(q), ' periods ------'], strvcat('Variable','        Correlation'), NamesNBCs, NBCs_NIPA, 10,12,3);
disp('    ');
dyntable(['------ International Correlations, news are ', num2str(q), ' periods ------'], strvcat('Variable','             Correlation'), NamesIBCs, IBCs_NIPA, 10,12,3);
disp('    ');
dyntable(['------ Correlation of RER with..., news are ', num2str(q), ' periods ------'], strvcat('Variable','             Correlation'), NamesRERcorrs, RERcorrs_NIPA, 10,12,3);
disp('    ');
*/