% ---------------------------- Weather Disasters -- Un-normalized -------------
application_str = 'Disasters';
application_dir = '/Users/mwatson/Dropbox/TVExtreme/ReplicationFiles/Stan/WeatherDamages/';
stan_draws_dir = [application_dir 'CSV/'];
model_str = 'Exceedance_GEV_4RW_Model';
data_fname = [matdir 'Disasters.mat'];
xlim_dates = [1975 2025];  % For Plotting
% Read in Data to get calendar dates
load(data_fname);
T = size(calvec,1);
tau = ones(T,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);
T = size(xi_draws,2);

% 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);
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

xi_pct_dam = xi_pct;
GEV_90_pct_dam = GEV_90_pct;
GEV_90_mean_dam = GEV_90_mean;

% ---------------------------- Weather Disasters -- Normalized -------------
application_str = 'Disasters_Normalized';
application_dir = '/Users/mwatson/Dropbox/TVExtreme/ReplicationFiles/Stan/WeatherDamages/';
stan_draws_dir = [application_dir 'CSV/'];
model_str = 'Exceedance_GEV_4RW_Model';
data_fname = [matdir 'Disasters_norm.mat'];
xlim_dates = [1975 2025];  % For Plotting
% Read in Data to get calendar dates
load(data_fname);
T = size(calvec,1);
tau = tau_norm;

% ------------------------------------------------------------

% 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);
T = size(xi_draws,2);

% 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_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

xi_pct_norm = xi_pct;
GEV_90_pct_norm = GEV_90_pct;
GEV_90_mean_norm = GEV_90_mean;
