% Summarize Results for 4RW model for City Size and Firm Size

% -------- City Size -------- 
application_str = 'CitySize';
application_dir = '/Users/mwatson/Dropbox/TVExtreme/ReplicationFiles/Stan/CitySize/';
stan_draws_dir = [application_dir 'CSV/'];
model_str = 'GEV_4RW_Model';

% Data
data_fname = [matdir 'CitySize_100.mat'];
load(data_fname);
T = size(cs_data,1);
calvec = yr;

tmp = readmatrix([stan_draws_dir application_str '.' model_str '.Draws.' 'xi' '.csv']);
xi_draws = tmp(2:end,2:end);  % Eliminate first row and first column 
tmp = readmatrix([stan_draws_dir application_str '.' model_str '.Draws.' 'sigma' '.csv']);
sigma_draws = tmp(2:end,2:end);  % Eliminate first row and first column 
tmp = readmatrix([stan_draws_dir application_str '.' model_str '.Draws.' 'mu' '.csv']);
mu_draws = tmp(2:end,2:end);  % Eliminate first row and first column 

nrep = size(xi_draws,1);
T1 = size(xi_draws,2);
if T1 ~= T
  error('T1 ~= T');
end

% GEV Draws 
GEV_90_draws = NaN(nrep,T);
for t = 1:T
    GEV_90_draws(:,t) = gevinv(0.90,xi_draws(:,t),sigma_draws(:,t),mu_draws(:,t));
end

% Compute Percentiles for time varying parameters
perc = [2.5 100*1/6 50 100*5/6 97.5];
n_perc = length(perc);
% GEV Parameters
mu_pct = NaN(T,n_perc);
sigma_pct = NaN(T,n_perc);
xi_pct = NaN(T,n_perc);
GEV_90_pct = NaN(T,n_perc);

for t = 1:T
  mu_pct(t,:) = prctile(mu_draws(:,t),perc);
  sigma_pct(t,:) = prctile(sigma_draws(:,t),perc);
  xi_pct(t,:) = prctile(xi_draws(:,t),perc);
  GEV_90_pct(t,:) = prctile(GEV_90_draws(:,t),perc);
end

% Save Key Results for plotting
xi_pct_cs = xi_pct;
GEV_90_pct_cs = GEV_90_pct;
yr_cs = yr;

% ---- Repeat for Firm Size ----
% application_str = 'EMP100';
application_str = 'FirmSize';
application_dir = '/Users/mwatson/Dropbox/TVExtreme/ReplicationFiles/Stan/FirmSize/';;
stan_draws_dir = [application_dir 'CSV/'];
model_str = 'GEV_4RW_Model';

% Data
str_data = [matdir 'FirmSize_top100_byYear.mat'];
load(str_data)
emp_all = emp_largest_normalized;

% Select years
yr = calvec;
yr_use = [1950 1973 1996 2019]';
yr = calvec;
emp_data = NaN(length(yr_use),size(emp_all,2));
for i = 1:length(yr_use)
    yr_idx = find(yr==yr_use(i));
    emp_data(i,:) = emp_all(yr_idx,:);
end
yr = yr_use;
T = size(emp_data,1);
calvec = yr;

tmp = readmatrix([stan_draws_dir application_str '.' model_str '.Draws.' 'xi' '.csv']);
xi_draws = tmp(2:end,2:end);  % Eliminate first row and first column 
tmp = readmatrix([stan_draws_dir application_str '.' model_str '.Draws.' 'sigma' '.csv']);
sigma_draws = tmp(2:end,2:end);  % Eliminate first row and first column 
tmp = readmatrix([stan_draws_dir application_str '.' model_str '.Draws.' 'mu' '.csv']);
mu_draws = tmp(2:end,2:end);  % Eliminate first row and first column 

nrep = size(xi_draws,1);
T1 = size(xi_draws,2);
if T1 ~= T
  error('T1 ~= T');
end

% GEV Draws 
GEV_90_draws = NaN(nrep,T);
for t = 1:T
    GEV_90_draws(:,t) = gevinv(0.90,xi_draws(:,t),sigma_draws(:,t),mu_draws(:,t));
end

% Compute Percentiles for time varying parameters
perc = [5 100*1/6 50 100*5/6 95];
n_perc = length(perc);
% GEV Parameters
mu_pct = NaN(T,n_perc);
sigma_pct = NaN(T,n_perc);
xi_pct = NaN(T,n_perc);
GEV_90_pct = NaN(T,n_perc);

for t = 1:T
  mu_pct(t,:) = prctile(mu_draws(:,t),perc);
  sigma_pct(t,:) = prctile(sigma_draws(:,t),perc);
  xi_pct(t,:) = prctile(xi_draws(:,t),perc);
  GEV_90_pct(t,:) = prctile(GEV_90_draws(:,t),perc);
end

% Save Key Results for plotting
xi_pct_fs = xi_pct;
GEV_90_pct_fs = GEV_90_pct;
yr_fs = yr;


