clear all;
small = 1.0e-10;
pinv_tol = 1.0e-05;
big = 1.0e+8;

% include path for utilities/functions
addpath('../Matlab_Functions');

% -- File Directories   
outdir = 'out/';
figdir = 'fig/';
matdir = 'mat/';

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

GEV_90_mean = NaN(T,1);

for t = 1:T
  GEV_90_mean(t) = mean(GEV_90_draws(:,t));
end

GEV_90_mean_unnorm = 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;

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

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

GEV_90_mean = NaN(T,1);

for t = 1:T
  GEV_90_mean(t) = mean(GEV_90_draws(:,t));
end

GEV_90_mean_norm = GEV_90_mean;

[calvec GEV_90_mean_unnorm GEV_90_mean_norm]
