%% Bailey, N., Holly, S. and Pesaran, H. P. (2015). A Two Stage Approach to Spatio-Temporal Analysis with Strong and Weak Cross-Sectional Dependence. 
                                                    %Journal of Applied Econometrics, forthcoming.
clear all

%% Input Data
m_series  = csvread('csv_residuals.csv') * 100; % Residuals from regressions using the cross-sectional averages approach
N = size(m_series, 1);
m_Cm = csvread('csv_Wmatrix_negative_connections.csv') - eye(N); % W_- (estimated using residuals from CSA regressions) 
m_Cp = csvread('csv_Wmatrix_positive_connections.csv') - eye(N); % W_+ (estimated using residuals from CSA regressions)

%% Row normalisation of W matrices
v_rowsum_Cm = sum(m_Cm, 2);
v_rowsum_Cp = sum(m_Cp, 2);
v_set = zeros(N, 1);
for i = 1:N
     if     (v_rowsum_Cm(i) ~= 0) && (v_rowsum_Cp(i) ~= 0); v_set(i) = 1;
     elseif (v_rowsum_Cm(i) ~= 0) && (v_rowsum_Cp(i) == 0); v_set(i) = 2;
     elseif (v_rowsum_Cm(i) == 0) && (v_rowsum_Cp(i) ~= 0); v_set(i) = 3;
     elseif (v_rowsum_Cm(i) == 0) && (v_rowsum_Cp(i) == 0); v_set(i) = 4;
     end
end
v_msa_without_neighm = find(v_rowsum_Cm == 0);
v_msa_without_neighp = find(v_rowsum_Cp == 0);
v_rowsum_Cm(v_msa_without_neighm) = 1;
v_rowsum_Cp(v_msa_without_neighp) = 1;

m_rowsum_Cm = repmat(v_rowsum_Cm, 1, N);
m_rowsum_Cp = repmat(v_rowsum_Cp, 1, N);
m_Wm = m_Cm ./ m_rowsum_Cm;
m_Wp = m_Cp ./ m_rowsum_Cp;

%% MLE estimation - based on Aquaro, M., Bailey, N. and Pesaran, M.H. (2014). Quasi-Maximum Likelihood Estimation of Spatial...
                                                                             %Models with Heterogeneous Coefficients. Under preparation
m_y = m_series(:, 2:end);
m_ly = m_series(:, 1:(end - 1));
m_lysm = m_Wm * m_ly;
m_lysp = m_Wp * m_ly;
a_x = ones([size(m_y) 4]);
a_x(:, :, 2) = m_ly;
a_x(:, :, 3) = m_lysm;
a_x(:, :, 4) = m_lysp;

% Choice of initial value
% 1- concentrated MLE: initial point for the maximization of the full
% likelihood function (reduces computational time)
v_initial_point = NaN;
concentrated = 1;
[m_eml m_vml] = fn_ml_N2psi_NKbeta_Nsgm_with_selection(m_y, a_x, m_Wm, m_Wp, concentrated, v_initial_point, v_set);
v_psim = m_eml(:, 1);
v_psip = m_eml(:, 2);
m_beta = m_eml(:, 3:(end - 1));
v_sgm  = m_eml(:, end);
m_beta_tr = m_beta';
v_beta = m_beta_tr(:);
v_initial_point = [v_psim; v_psip; v_beta; v_sgm];

% 2- full (non-concentrated) MLE
concentrated = 0;
[m_eml m_vml] = fn_ml_N2psi_NKbeta_Nsgm_with_selection(m_y, a_x, m_Wm, m_Wp, concentrated, v_initial_point, v_set);

%% Matrix of estimates per MSA and significance indicators
% Estimates rearranged in line with Tables 7 and 8 of BHP paper (lambda1, psi0+, psi0-, psi1+, psi1-, sigma_zeta)
m_eml = m_eml(:, [4 2 1 6 5 7]);
m_vml = m_vml(:, [4 2 1 6 5 7]);

% Significance: testing
m_set = Inf(size(m_eml));
for i_row = 1:size(m_eml)
     if     (v_set(i_row) == 1) m_set(i_row, :) = [1 1 1 1 1 1]; % note: gamma, psi0+, psi0-, psi1+, psi1-, sgm
     elseif (v_set(i_row) == 2) m_set(i_row, :) = [1 0 1 0 1 1];
     elseif (v_set(i_row) == 3) m_set(i_row, :) = [1 1 0 1 0 1];
     elseif (v_set(i_row) == 4) m_set(i_row, :) = [1 0 0 0 0 1];
     end
end

cval = 1.96; %cc%
m_reject = Inf(size(m_eml));
for i_row = 1:size(m_eml, 1)
     for i_col = 1:size(m_eml, 2)
          if m_set(i_row, i_col) == 1
               tval = m_eml(i_row, i_col) / sqrt(m_vml(i_row, i_col));
               m_reject(i_row, i_col) = abs(tval) > cval;
          elseif m_set(i_row, i_col) == 0
               m_reject(i_row, i_col) = -1;
          end
     end
end

%% Output of MLE estimates and significance indicators per MSA
disp('MLE estimates');
disp(num2str(m_eml));
disp(' ');
disp('MLE significance indicators');
disp(num2str(m_reject));
disp(' ');

%% Output of Tables 7 and 8 in Section 4.3 of BHP paper 
m_summ_stat = zeros(5, size(m_eml, 2));
for i_col = 1:size(m_eml, 2)
     v_est_temp = m_eml(:, i_col);
     v_rej_temp = m_reject(:, i_col);
     if ((i_col == 3) || (i_col == 5))
          v_est = v_est_temp(find((v_set == 1) | (v_set == 2)), 1);
          v_rej = v_rej_temp(find((v_set == 1) | (v_set == 2)), 1);
     elseif ((i_col == 2) || (i_col == 4))
          v_est = v_est_temp(find((v_set == 1) | (v_set == 3)), 1);
          v_rej = v_rej_temp(find((v_set == 1) | (v_set == 3)), 1);
     elseif ((i_col == 1) || (i_col == 6))
          v_est = v_est_temp;
          v_rej = v_rej_temp;
     end
     m_summ_stat(1, i_col) = quantile(v_est, .5);
     m_summ_stat(2, i_col) = mean(v_est);
     m_summ_stat(3, i_col) = std(v_est) / sqrt(length(v_est));
     m_summ_stat(4, i_col) = mean(v_rej) * 100;
     m_summ_stat(5, i_col) = length(v_est);
end

texfilename = 'in_tb7.tex';
zz = fopen(texfilename, 'w');
fprintf(zz, '%% -------- Table 7 -------- %%\n');
mat = m_summ_stat;
for i_row = 1:size(mat, 1)
     if i_row == 1
          fprintf(zz, 'Median                        & ');
          for i_col = 1:size(mat, 2)
               if i_col < size(mat, 2)
                    fprintf(zz, '%8.4f & ', mat(i_row, i_col));
               else
                    fprintf(zz, '%8.4f \\\\ \n', mat(i_row, i_col));
               end
          end
     elseif i_row == 2
          fprintf(zz, 'Mean Group Estimates          & ');
          for i_col = 1:size(mat, 2)
               if i_col < size(mat, 2)
                    fprintf(zz, '%8.4f & ', mat(i_row, i_col));
               else
                    fprintf(zz, '%8.4f \\\\ \n', mat(i_row, i_col));
               end
          end
     elseif i_row == 3
          fprintf(zz, '                              & ');
     for i_col = 1:size(mat, 2)
          if i_col < size(mat, 2)
               fprintf(zz, '(%6.4f) & ', mat(i_row, i_col));
          else
               fprintf(zz, '(%6.4f) \\\\ \n', mat(i_row, i_col));
          end
     end
     elseif i_row == 4
          fprintf(zz, '\\%% significant (at 5\\%% level) & ');
          for i_col = 1:size(mat, 2)
               if i_col < size(mat, 2)
                    fprintf(zz, '%6.1f\\%% & ', mat(i_row, i_col));
               else
                    fprintf(zz, '       - \\\\ \n');
               end
          end
     elseif i_row == 5
          fprintf(zz, 'Number of non-zero coef.      & ');
          for i_col = 1:size(mat, 2)
               if i_col < size(mat, 2)
                    fprintf(zz, '%8.0f & ', mat(i_row, i_col));
               else
                    fprintf(zz, '%8.0f \\\\ \n', mat(i_row, i_col));
               end
          end
     end
end
fclose(zz);
fprintf('\nTable 7. Output saved in the tex file: "%s"\n', texfilename)

texfilename = 'in_tb8.tex';
zz = fopen(texfilename, 'w');
fprintf(zz, '%% -------- Table 8 -------- %%\n');
qwe = importdata('cvs_dictionary.csv', ',');
qwe.textdata(1, :) = []; % kill header
v_region = qwe.textdata(:, 1);
v_region_names = unique(v_region);
clear qwe

for i_reg = [4 3 7 2 5 8 6 1]
     region_name = v_region_names(i_reg);
     v_ind = find(ismember(v_region, region_name) == 1);
     m_eml_temp = m_eml(v_ind, :);
     m_reject_temp = m_reject(v_ind, :);
     v_set_temp = v_set(v_ind, :);

     m_summ_stat = zeros(5, size(m_eml_temp, 2));
     for i_col = 1:size(m_eml_temp, 2)
          v_est_temp = m_eml_temp(:, i_col);
          v_rej_temp = m_reject_temp(:, i_col);
          if ((i_col == 3) || (i_col == 5))
               v_est = v_est_temp(find((v_set_temp == 1) | (v_set_temp == 2)), 1);
               v_rej = v_rej_temp(find((v_set_temp == 1) | (v_set_temp == 2)), 1);
          elseif ((i_col == 2) || (i_col == 4))
               v_est = v_est_temp(find((v_set_temp == 1) | (v_set_temp == 3)), 1);
               v_rej = v_rej_temp(find((v_set_temp == 1) | (v_set_temp == 3)), 1);
          elseif ((i_col == 1) || (i_col == 6))
               v_est = v_est_temp;
               v_rej = v_rej_temp;
          end
          m_summ_stat(1, i_col) = quantile(v_est, .5);
          m_summ_stat(2, i_col) = mean(v_est);
          m_summ_stat(3, i_col) = std(v_est) / sqrt(length(v_est));
          m_summ_stat(4, i_col) = mean(v_rej) * 100;
          m_summ_stat(5, i_col) = length(v_est);
     end
     fprintf(zz, '\n%%%s\n', char(region_name));
     mat = m_summ_stat;
     for i_row = 1:size(mat, 1)
          if i_row == 1
               fprintf(zz, 'Median                        & ');
               for i_col = 1:size(mat, 2)
                    if i_col < size(mat, 2)
                         fprintf(zz, '%8.4f & ', mat(i_row, i_col));
                    else
                         fprintf(zz, '%8.4f \\\\ \n', mat(i_row, i_col));
                    end
               end
          elseif i_row == 2
               fprintf(zz, 'Mean Group Estimates          & ');
               for i_col = 1:size(mat, 2)
                    if i_col < size(mat, 2)
                         fprintf(zz, '%8.4f & ', mat(i_row, i_col));
                    else
                         fprintf(zz, '%8.4f \\\\ \n', mat(i_row, i_col));
                    end
               end
          elseif i_row == 3
               fprintf(zz, '                              & ');
          for i_col = 1:size(mat, 2)
               if i_col < size(mat, 2)
                    fprintf(zz, '(%6.4f) & ', mat(i_row, i_col));
               else
                    fprintf(zz, '(%6.4f) \\\\ \n', mat(i_row, i_col));
               end
          end
          elseif i_row == 4
               fprintf(zz, '\\%% significant (at 5\\%% level) & ');
               for i_col = 1:size(mat, 2)
                    if i_col < size(mat, 2)
                         fprintf(zz, '%6.1f\\%% & ', mat(i_row, i_col));
                    else
                         fprintf(zz, '       - \\\\ \n');
                    end
               end
          elseif i_row == 5
               fprintf(zz, 'Number of non-zero coef.      & ');
               for i_col = 1:size(mat, 2)
                    if i_col < size(mat, 2)
                         fprintf(zz, '%8.0f & ', mat(i_row, i_col));
                    else
                         fprintf(zz, '%8.0f \\\\ \n', mat(i_row, i_col));
                    end
               end
          end
     end
end
fclose(zz);
fprintf('\nTable 8. Output saved in the tex file: "%s"\n', texfilename)
