% Forecast Results

% Descriptive Statists
clear all;
small = 1.0e-10;

% Model Name
mod_name = '_baseline';
% mod_name = '_fgamma05';
% mod_name = '_fgamma15';
% mod_name = '_q_8';
% mod_name = '_q_20';

% -- File Directories   
figdir = 'fig/';
outdir = 'out/';
matdir = 'mat/';
mcmcdir = 'mcmc/';
figdir_fcst = 'fig_fcst/';

fstr = [matdir 'country_code']; load(fstr);
fstr = [matdir 'country']; load(fstr);
n_c = size(country_code,1);

fstr = [matdir 'yfit_mat'];load(fstr);  %Trends .. for plotting

% load Paths
fstr = [mcmcdir 'path_F_draws' mod_name];load(fstr);
fstr = [mcmcdir 'paths_U_draws' mod_name];load(fstr);
fstr = [matdir 'yp']; load(fstr);
fstr = [matdir 'pop']; load(fstr);

% Rename;
u_draws = paths_U_draws;
f_draws = path_F_draws;
clear paths_U_draws;
clear path_F_draws;

T = 118;  % Number of insample obs

% Compute Growth Rate Paths
U_gr_paths_50 = squeeze(u_draws(T+50,:,:)-u_draws(T,:,:))/50;
U_gr_paths_100 = squeeze(u_draws(T+100,:,:)-u_draws(T,:,:))/100;
f_gr_paths_50 = (f_draws(T+50,:)-f_draws(T,:))/50;
f_gr_paths_100 = (f_draws(T+100,:)-f_draws(T,:))/100;
Y_gr_paths_50 = U_gr_paths_50 + repmat(f_gr_paths_50,n_c,1);
Y_gr_paths_100 = U_gr_paths_100 + repmat(f_gr_paths_100,n_c,1);

% Construct Growth Rate Prediction intervals for each country
% Save Results 
outfile_name = [outdir 'Table_5' mod_name '.csv'];
fileID = fopen(outfile_name,'w');
fprintf(fileID,'Forecast Growth Rates (PAAR) \n');
fprintf(fileID,',,50Years,,,,,,,100Years \n');
pct = [0.05 1/6 0.50 5/6 0.95]';
str1 = ['Mean,.05,.16,.50,.84,.95'];
fprintf(fileID,[',,' str1 ',,' str1 '\n']);
fprintf(fileID,'F,,');
tmp = pctile_mean(f_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(f_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f','\n');

% Use population weights for growth rates
fprintf(fileID,'All,,');
pop_g = pop(end,:)';     % All countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50;
g_gr_paths_100 = wght'*Y_gr_paths_100;
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

% Load Lists
load([matdir 'oecd_list']);
i_oecd = zeros(n_c,1);
nn = size(oecd_list,1);
for j = 1:nn;
    ss = char(oecd_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_oecd(jj) = 1;
    end;
end;

ii = i_oecd==1;
fprintf(fileID,'\n OECD,,');
pop_g = pop(end,ii==1)';     % OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

ii = i_oecd==0;
fprintf(fileID,'\n Non-OECD,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

ilist_Asia1 = {'THA';'SGP';'TWN';'HKG';'KOR';'MYS'};
ilist_Asia2 = {'MLI';'MMR';'LKA';'CHN';'IDN';'VNM';'LAO';'IND';'PAN'};
ilist_Soviet = {'HUN';'SRB';'BGR';'HRV';'RUS';'ROU'};
ilist_AS = {'AUS';'CAN';'GBR';'NZL';'USA'};
ilist_Europe = {'FIN';'CHE';'NOR';'DNK';'SWE';'NLD';'ITA';'BEL';'FRA';'GRC'};
ilist_last = {'DEU';'AUT';'ESP';'PRT';'SYR';'YEM';'POL';'DOM';'BIH';'CHL';'COG';'DZA';'EGY';'GTM';'IRL';'IRN';'IRQ';'ISR';'JOR';'JPN';'LBN';'LBR';'MEX';'MNG';'OMN';'PHL';'PRY';'PSE';'SAU';'TUN';'TUR';'ZAF'};

g_list = ilist_Asia1;
i_g = zeros(n_c,1);
nn = size(g_list,1);
for j = 1:nn;
    ss = char(g_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_g(jj) = 1;
    end;
end;
i_Asia1 = i_g;

g_list = ilist_Asia2;
i_g = zeros(n_c,1);
nn = size(g_list,1);
for j = 1:nn;
    ss = char(g_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_g(jj) = 1;
    end;
end;
i_Asia2 = i_g;

g_list = ilist_AS;
i_g = zeros(n_c,1);
nn = size(g_list,1);
for j = 1:nn;
    ss = char(g_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_g(jj) = 1;
    end;
end;
i_AS = i_g;

g_list = ilist_Europe;
i_g = zeros(n_c,1);
nn = size(g_list,1);
for j = 1:nn;
    ss = char(g_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_g(jj) = 1;
    end;
end;
i_Europe = i_g;

g_list = ilist_Soviet;
i_g = zeros(n_c,1);
nn = size(g_list,1);
for j = 1:nn;
    ss = char(g_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_g(jj) = 1;
    end;
end;
i_Soviet = i_g;

g_list = ilist_last;
i_g = zeros(n_c,1);
nn = size(g_list,1);
for j = 1:nn;
    ss = char(g_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_g(jj) = 1;
    end;
end;
i_last = i_g;

ii = i_Asia1+i_Asia2+i_AS+i_Europe+i_Soviet+i_last;
i_other = ii==0;

ii = i_Asia1==1;
sum(ii)
fprintf(fileID,'\n Asia1,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
 
ii = i_Asia2==1;
sum(ii)
fprintf(fileID,'\n Asia2,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

ii = i_AS==1;
sum(ii)
fprintf(fileID,'\n AS,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

ii = i_Europe==1;
sum(ii)
fprintf(fileID,'\n Europe,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

ii = i_Soviet==1;
sum(ii)
fprintf(fileID,'\n Soviet,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

ii = i_other==1;
sum(ii)
fprintf(fileID,'\n Other,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

ii = i_last==1;
sum(ii)
fprintf(fileID,'\n last,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

% Compute and Save Country-specific forecasts
% Construct Growth Rate Prediction intervals for each country
% Save Results 
outfile_name = [outdir 'Table_10_detail_byCountry' mod_name '.csv'];
fileID = fopen(outfile_name,'w');
fprintf(fileID,'Forecast Growth Rates (PAAR) \n');
fprintf(fileID,',,50Years,,,,,,,100Years \n');
pct = [0.05 1/6 0.50 5/6 0.95]';
str1 = ['Mean,.05,.16,.50,.84,.95'];
fprintf(fileID,[',,2017Value,,' str1 ',,' str1 '\n']);
fprintf(fileID,'F,,');
tmp = mean(f_draws(T,:)');
fprintf(fileID,'%5.2f ,,',tmp);
tmp = pctile_mean(f_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(f_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f','\n');

fprintf(fileID,'\n\n');
for ic = 1:n_c;
    str = char(country_code(ic));
    fprintf(fileID,[str ',,']);
    tmp1 = squeeze(u_draws(T,ic,:));
    tmp = mean(f_draws(T,:)') + mean(tmp1');
    fprintf(fileID,'%5.2f ,,',tmp);
    tmp = pctile_mean(Y_gr_paths_50(ic,:)',pct)';
    prtmat_comma(100*tmp,fileID,'%5.2f',',,');
    tmp = pctile_mean(Y_gr_paths_100(ic,:)',pct)';
    prtmat_comma(100*tmp,fileID,'%5.2f',',,');
    
    % Construct Median Log-level after 100 years
    ytmp = squeeze(u_draws(end,ic,:))+f_draws(end,:)';
    med_100 = median(ytmp);
    fprintf(fileID,'%5.2f \n',med_100);
end;

% Save Some Results -- used below
post_m_perc_save_multivariate = NaN(6,2,n_c);
for ic = 1:n_c;
    post_m_perc_save_multivariate(:,1,ic) = pctile_mean(Y_gr_paths_50(ic,:)',pct);
    post_m_perc_save_multivariate(:,2,ic) = pctile_mean(Y_gr_paths_100(ic,:)',pct);
end;
post_f_perc_save = NaN(6,2);
post_f_perc_save(:,1) = pctile_mean(f_gr_paths_50(1,:)',pct);
post_f_perc_save(:,2) = pctile_mean(f_gr_paths_100(1,:)',pct);


% Figure 12 -- Multivariate Prediction intervals by country
M50 = squeeze(post_m_perc_save_multivariate(2:6,1,:));
M100 = squeeze(post_m_perc_save_multivariate(2:6,2,:));

% Sort by average YP over last 10 years
yp_avg = mean(yp(end-4:end,:))';
[~,II] = sort(yp_avg);
country_code_srt = country_code(II);
M50 = M50(:,II)';
M100 = M100(:,II)';

% Error bands for plotting
trend = (1:1:n_c)';
trendp = trend+0.1;


% List country names for future use;
outfile_name = [outdir 'plot_univariat_multivariate_fcst_intervals_country_names' mod_name '.csv'];
fileID = fopen(outfile_name,'w');
for i = 1:n_c;
    fprintf(fileID,'%3i ,',i);
    str = char(country_code_srt(i));
    fprintf(fileID,[str '\n']);
end;

figure_12;

% Figure 13
% Read in Univariate results are compute predictive percentiles
U50_pct = NaN(n_c,5);
U100_pct = NaN(n_c,5);
% Read in UMs values
umdir = ['/Users/mwatson/Dropbox/Shared_Folders/SCC/output_matlab/ndata/'];
fstr = [umdir 'single_paths_baseline']; load(fstr);

% Compute univariate forecast intervals and save them
for ic = 1:n_c;
    pc = squeeze(paths(:,ic,:));
    g50 = 100*(pc(T+50,:)-pc(T,:))/50;
    tmp = pctile(g50',pct)';
    U50_pct(ic,:) = tmp;
    g100 = 100*(pc(T+100,:)-pc(T,:))/100;
    tmp = pctile(g100',pct)';
    U100_pct(ic,:) = tmp;
end;

% Save Univariate percentiles to a file
outfile_name = [outdir 'Univariate_Forecast_Intervals.csv'];
fileID = fopen(outfile_name,'w');
fprintf(fileID,'Country 50 the 100 year percentiles \n');
for ic = 1:n_c;
    str = char(country_code(ic));
    fprintf(fileID,[str ',,']);
    prtmat_comma(U50_pct(ic,:),fileID,'%5.3f',',,');
    prtmat_comma(U100_pct(ic,:),fileID,'%5.3f',',,');
    
    % Construct Median Log-level after 100 years
    ytmp = squeeze(paths(end,ic,:));
    med_100 = median(ytmp);
    fprintf(fileID,'%5.2f \n',med_100);
    
end;

M50 = squeeze(post_m_perc_save_multivariate(2:6,1,:));
M100 = squeeze(post_m_perc_save_multivariate(2:6,2,:));

% Sort by average YP over last 10 years
yp_avg = mean(yp(end-4:end,:))';
[~,II] = sort(yp_avg);
country_code_srt = country_code(II);
M50 = M50(:,II)';
M100 = M100(:,II)';
U50 = U50_pct(II,:)/100;
U100 = U100_pct(II,:)/100;

% Try error bands for plotting
trend = (1:1:n_c)';
trendp = trend+0.25;

figure_13;

%   Figure 11
% Plot Individual figures for each country (save for appendix)
ffcst_pct = NaN(T+100,5);
for t = 1:T+100;
    ffcst_pct(t,:) = pctile(f_draws(t,:)',pct');
end;
yfcst_m_pct = NaN(T+100,n_c,5);
yfcst_u_pct = NaN(T+100,n_c,5);
for ic = 1:n_c;
    U_path = squeeze(u_draws(:,ic,:));
    Y_path = U_path + f_draws;
    Y_u_path = squeeze(paths(:,ic,:));
    for t = 1:T+100;
      yfcst_m_pct(t,ic,:) = pctile(Y_path(t,:)',pct');
      yfcst_u_pct(t,ic,:) = pctile(Y_u_path(t,:)',pct');
    end;
end;

calvec_p = (1900:2117)';
yall = NaN(T+100,n_c);
yall(1:T,:) = yfit_mat(:,1:113);
fall = NaN(T+100,1);
fall(1:T) = ffcst_pct(1:T,3);
fcst = NaN(T+100,5);
fcst(T+1:end,:) = ffcst_pct(T+1:end,:);

% Select Countries for 9 panel plot
%country_list = {'USA';'JPN';'SGP';'ESP';'HUN';'BRA';'CHN';'LBR'};
country_list = {'USA';'JPN';'SGP';'ESP';'HUN';'CHN';'SYR';'LBR'};

figure_11;


% Individual Panel Plots (one for each country)
figure_11_individual_panels;

% Figure 14
% Select Countries for 6 panel plot
country_list = {'CAF';'DNK';'GRC';'IND';'KOR';'ECU'};
%country_list = {'CAF';'DNK';'GRC';'IND';'YEM';'ECU'};
figure_14;

% Individual Panel Plots (one for each country)
figure_14_individual_panels;
