% Summarize Results for 4RW model

% ---------------------------- VW Returns -- Largest -------------
application_str = 'Returns_Large';
application_dir = '/Users/mwatson/Dropbox/TVExtreme/ReplicationFiles/Stan/Returns/';
stan_draws_dir = [application_dir 'CSV/'];
model_str = 'GEV_4RW_Model';
data_fname =  [matdir 'Returns_Large.mat'];
% Read in Data to get calendar dates
load(data_fname);
T = length(ret);
xlim_dates = [1920 2025];  % For Plotting
sign_scale_plot = 1.0;   % 1.0 for positive returns, -1.0 for negative returns
% % ------------------------------------------------------------

% Load Draws
% Constant Parameters
tmp = readmatrix([stan_draws_dir application_str '.' model_str '.Draws.' 'gamma_alpha' '.csv']); gamma_alpha_draws = tmp(2:end,2:end);  % Eliminate first row and first column 
tmp = readmatrix([stan_draws_dir application_str '.' model_str '.Draws.' 'gamma_xi' '.csv']); gamma_xi_draws = tmp(2:end,2:end);  % Eliminate first row and first column 
tmp = readmatrix([stan_draws_dir application_str '.' model_str '.Draws.' 'gamma_s' '.csv']); gamma_s_draws = tmp(2:end,2:end);  % Eliminate first row and first column 
tmp = readmatrix([stan_draws_dir application_str '.' model_str '.Draws.' 'gamma_m' '.csv']); gamma_m_draws = tmp(2:end,2:end);  % Eliminate first row and first column 

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

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);
GEV_90_mean = NaN(T,1);

parfor 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);
  GEV_90_mean(t) = mean(GEV_90_draws(:,t));
end

% Scale mu and GEV for plotting
mu_pct = sign_scale_plot*mu_pct;
GEV_90_pct = sign_scale_plot*GEV_90_pct;
GEV_90_mean = sign_scale_plot*GEV_90_mean;

xi_pct_large = xi_pct;
GEV_90_pct_large = GEV_90_pct;
GEV_90_mean_large = GEV_90_mean;


% % ---------------------------- VW Returns -- Smallest -------------
application_str = 'Returns_Small';
pplication_dir = '/Users/mwatson/Dropbox/TVExtreme/ReplicationFiles/Stan/Returns/';
stan_draws_dir = [application_dir 'CSV/'];
model_str = 'GEV_4RW_Model';
data_fname =  [matdir 'Returns_Small.mat'];
% Read in Data to get calendar dates
load(data_fname);
T = length(ret);
xlim_dates = [1920 2025];  % For Plotting
sign_scale_plot = -1.0;   % 1.0 for positive returns, -1.0 for negative returns
marker_size = 1;
% ------------------------------------------------------------

% Load Draws
% Constant Parameters
tmp = readmatrix([stan_draws_dir application_str '.' model_str '.Draws.' 'gamma_alpha' '.csv']); gamma_alpha_draws = tmp(2:end,2:end);  % Eliminate first row and first column 
tmp = readmatrix([stan_draws_dir application_str '.' model_str '.Draws.' 'gamma_xi' '.csv']); gamma_xi_draws = tmp(2:end,2:end);  % Eliminate first row and first column 
tmp = readmatrix([stan_draws_dir application_str '.' model_str '.Draws.' 'gamma_s' '.csv']); gamma_s_draws = tmp(2:end,2:end);  % Eliminate first row and first column 
tmp = readmatrix([stan_draws_dir application_str '.' model_str '.Draws.' 'gamma_m' '.csv']); gamma_m_draws = tmp(2:end,2:end);  % Eliminate first row and first column 

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

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);
GEV_90_mean = NaN(T,1);

parfor 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);
  GEV_90_mean(t) = mean(GEV_90_draws(:,t));
end

% Scale mu and GEV for plotting
mu_pct = sign_scale_plot*mu_pct;
GEV_90_pct = sign_scale_plot*GEV_90_pct;
GEV_90_mean = sign_scale_plot*GEV_90_mean;

xi_pct_small = xi_pct;
GEV_90_pct_small = GEV_90_pct;
GEV_90_mean_small = GEV_90_mean;





