% BS paper: model non-linearised, letting dynare do it - simple linearisation, not log-linearisation

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

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

var p rerCONS p2T c_rel nx1_y1 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}

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;
pi1 = (c1 - PSI_1*n1^ni_1)^(gamma_1-1) ;
pi2 = (c2 - PSI_2*n2^ni_2)^(gamma_2-1) ;

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

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

% new variable for investment adjustment cost function: sI
%sI1 = kappa_1*(i1/i1(-1) - 1)^2/2;
%sI2 = kappa_2*(i2/i2(-1) - 1)^2/2;
%sprimeI1 = kappa_1*(i1/i1(-1) - 1);
%sprimeI2 = kappa_2*(i2/i2(-1) - 1);

phi1 = lamda1*(1 - kappa_1*(i1/i1(-1) - 1)^2/2 - i1*kappa_1*(i1/i1(-1) - 1)/i1(-1)) + beta_1*lamda1(+1)*kappa_1*(i1(+1)/i1 - 1)*(i1(+1)/i1)^2;
phi2 = lamda2*(1 - kappa_2*(i2/i2(-1) - 1)^2/2 - i2*kappa_2*(i2/i2(-1) - 1)/i2(-1)) + beta_2*lamda2(+1)*kappa_2*(i2(+1)/i2 - 1)*(i2(+1)/i2)^2;
zita1 = lamda1*deltaprimeU1 ;
zita2 = lamda2*deltaprimeU2 ;

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

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

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


% Stochastic processes 
% if you express it like this, the initial condition on Z=0 and the production function should be expressed
% as exp(zT1)*ST1^alphaT_1*nT1^(1-alphaT_1)
%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}) ;

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

% if you express it like this, the initial condition on Z=0 and the production function should be expressed
% as zT1*ST1^alphaT_1*nT1^(1-alphaT_1)
%log(zT1) = A1T_11*log(zT1(-1)) + A1T_12*log(zT2(-1)) + A1T_13*log(zN1(-1)) + A1T_14*log(zN2(-1)) + eT1(-@{q}) ;
%log(zT2) = A2T_11*log(zT1(-1)) + A2T_12*log(zT2(-1)) + A2T_13*log(zN1(-1)) + A2T_14*log(zN2(-1)) + eT2(-@{q}) ;
%log(zN1) = A1N_11*log(zT1(-1)) + A1N_12*log(zT2(-1)) + A1N_13*log(zN1(-1)) + A1N_14*log(zN2(-1)) + eN1(-@{q}) ;
%log(zN2) = A2N_11*log(zT1(-1)) + A2N_12*log(zT2(-1)) + A2N_13*log(zN1(-1)) + A2N_14*log(zN2(-1)) + eN2(-@{q}) ;

% Both ways of expressing the stochastic process/ production function produce the SAME results

% 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
mi1 = phi1*omega_1*a1^(-1/theta_1)*yT1^(1/theta_1);
mi1 = phi2*(1-omega_2)*a2^(-1/theta_2)*yT2^(1/theta_2);
mi2 = phi1*(1-omega_2)*b1^(-1/theta_1)*yT1^(1/theta_1);
mi2 = phi2*omega_2*b2^(-1/theta_2)*yT2^(1/theta_2);

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

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

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

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

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

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

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

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

% Investment
k1 = (1-deltaU1)*k1(-1) + i1*(1 - kappa_1*(i1/i1(-1) - 1)^2/2);
k2 = (1-deltaU2)*k2(-1) + i2*(1 - kappa_2*(i2/i2(-1) - 1)^2/2);

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

% Labour
n1 = nT1 + nN1 ;
n2 = nT2 + nN2 ;

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

% Equilibrium

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

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

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

% Other - definitions

% relative consumption
c_rel = c1/c2 ;
cN1_cT1 = cN1/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
y_rel = y1_real/y2_real;
% relative labour
n_rel = n1/n2 ;

mi_phi = mi1/phi1 ;
end;

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

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

rerNT = p2N_ss/p1N_ss;
rerLM = 1;
cN1_cT1 = cN1_ss/cT1_ss;
n_rel = n_ss_1/n_ss_1;
mi_phi=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 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;

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)));
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% This is a bit different than the one I have, cos I do not have a variable nx1_y1_real1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%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: 8 
Names = ['Consumption                 '; % position in var-cov matrix:  3
         'Investment                  '; % position in var-cov matrix: 12
         'Hours                       '; % position in var-cov matrix: 14
         '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:  20
         'NX over GDP (current prices)'];% position in var-cov matrix:  7


RelSDs = [sqrt(oo_.var(3,3)/oo_.var(8,8)) ; sqrt(oo_.var(12,12)/oo_.var(8,8)) ; sqrt(oo_.var(14,14)/oo_.var(8,8)) ; sqrt(oo_.var(2,2)/oo_.var(8,8)) ; sqrt(oo_.var(1,1)/oo_.var(8,8)) ; sqrt(oo_.var(20,20)/oo_.var(8,8)) ; sqrt(oo_.var(7,7)/oo_.var(8,8))];
% Within country co-movements
NBCs = [oo_.var(3,8)/(sqrt(oo_.var(3,3))*sqrt(oo_.var(8,8))) ; oo_.var(15,8)/(sqrt(oo_.var(15,15))*sqrt(oo_.var(8,8))) ; oo_.var(7,8)/(sqrt(oo_.var(7,7))*sqrt(oo_.var(8,8))) ];
NamesNBCs = ['Consumption                 '; 
             'Hours                       '; 
             'Investment                  '; 
             'NX over GDP (real)          '];

IBCs = [oo_.var(8,9)/(sqrt(oo_.var(8,8))*sqrt(oo_.var(9,9))) ; 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,7)/(sqrt(oo_.var(2,2))*sqrt(oo_.var(7,7))) ]; 
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(3)-(-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(3)-(-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(['irfsNL_', num2str(@{q}), ' = oo_.irfs;']); 
s = ['save irfsNL_' num2str(@{q}) '.mat, irfsNL_' num2str(@{q});] ;
eval(s) ;
%close all