% MLEs and Test Statistics for Disasters Data .. note needs to be run twice .. once for Raw data and then for Normalized data
clear all;
small = 1.0e-10;
pinv_tol = 1.0e-05;
big = 1.0e+8;

% Add new directory to path to find functions
addpath('../Matlab_Functions/');

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

% Set options for optimizers
options_fminunc = optimset('Display','off');
options_fminsearch = optimset('Display','off');

% Identifiers for files 
% -------------------------- Disaster Data: Raw --------------
% application_str = 'WeatherDamages';
% data_fname = [matdir 'Disasters.mat'];
% load(data_fname);
% % Rename some variables
% data_mat = y';
% nobs_vec = nobs;
% % --- Actual Data ---
% Y_data = data_mat;
% Y_nonmissing = Y_data(isnan(Y_data)==0);
% Y_nonmissing = Y_nonmissing(:);
% Y_data(Y_data==-1) = NaN;
% T = size(Y_data,1);
% % Value of threshold (tau) ... as a vector
% tau_vec = ones(T,1);
% % String for loading CVS adjustment parameters
% str_cvs = 'tau';
% % -------------------------------------------------------
% -------------------------- Disaster Data: Normalized --------------
application_str = 'WeatherDamages_Normalized';
data_fname = [matdir 'Disasters_norm.mat'];
load(data_fname);
% Rename some variables
data_mat = y_norm';
nobs_vec = nobs;
% --- Actual Data ---
Y_data = data_mat;
Y_nonmissing = Y_data(isnan(Y_data)==0);
Y_nonmissing = Y_nonmissing(:);
Y_data(Y_data==-1) = NaN;
T = size(Y_data,1);
% Value of threshold (tau) ... as a vector
tau_vec = tau_norm;
% String for loading CVS adjustment parameters
str_cvs = 'tvt';
% % -------------------------------------------------------

% Note -- parameterization is in terms of GEV parameters xsi, sigma and mu
xsi_start = 0.8;
sigma_start = 1.0;
mu_start = 0.4;

% Find MLE of GEV parameters using P-GP likelihood
x_start = [xsi_start sigma_start mu_start];
f_unc = @(x) minus_llf_exceedance_gev(x,Y_data,nobs_vec,tau_vec);
[x_1,fval_1]=fminunc(f_unc,x_start,options_fminunc); 
[x_2,fval_2] = fminsearch(f_unc,x_start,options_fminsearch); 
if fval_1 < fval_2
    x_mle = x_1;
    fval = fval_1;
else
    x_mle = x_2;
    fval = fval_2;
end
llf_unconstrained = -fval;
xsi_mle = x_mle(1);
sigma_mle = x_mle(2);
mu_mle = x_mle(3);

% Get Score and Information matrix
% % Compute Score and Hessian
[score_mat_unc,info_op_unc,info_h_unc] = score_info_Exceedance_GEV_tvtau(Y_data,nobs_vec,tau_vec,xsi_mle,sigma_mle,mu_mle);
score_mat_unc = score_mat_unc - mean(score_mat_unc);    % score_mat_unc should have mean zero ... impose this

info = info_h_unc;  % Use Hessian
cov_mle = inv(info)/T;
se_mle = sqrt(diag(cov_mle));
se_xsi_mle = se_mle(1);
se_sigma_mle = se_mle(2);
se_mu_mle = se_mle(3);

% Compute Nyblom Statistics -- All parameters
L_Stat = compute_nyblom(score_mat_unc,info_h_unc);
pv_L_Stat = pvalue_LStat(L_Stat,length(x_mle));
% Read in parameters for adjustment factor
fname = [cv_dir 'cvx_test0_' str_cvs '_' num2str(T) '.csv'];
c_parm = readmatrix(fname);
adj = exp(c_parm(1) + c_parm(2)*xsi_mle + c_parm(3)*xsi_mle^2);
adj_L_Stat = adj*L_Stat;


% Construct Test that mu = sigma/xsi
x_start = [xsi_mle sigma_mle];
f_con_mu = @(x) minus_llf_exceedance_gev_con_mu(x,Y_data,nobs_vec,tau_vec);
[x_1,fval_1] = fminunc(f_con_mu,x_start,options_fminunc);
[x_2,fval_2] = fminsearch(f_con_mu,x_start,options_fminsearch); 
if fval_1 < fval_2
  x_mle_constained_mu = x_1;
  fval = fval_1;
else
  x_mle_constained_mu = x_2;
  fval = fval_2;
end
llf_constrained_mu = -fval;
%  Compute adnjustment factor for test statistic
fname = [cv_dir 'cvx_test3_' str_cvs '_' num2str(T) '.csv'];
c_parm = readmatrix(fname);
adj = exp(c_parm(1) + c_parm(2)*xsi_mle + c_parm(3)*xsi_mle^2);
llf_dif = llf_unconstrained-llf_constrained_mu;
Stat_con_mu = adj*llf_dif; 

% Construct Test that mu = sigma/xsi and xsi = 1
x_start = [sigma_mle];
f_con_xsi_mu = @(x) minus_llf_exceedance_gev_con_xsi_mu(x,Y_data,nobs_vec,tau_vec);
[x_1,fval_1] = fminunc(f_con_xsi_mu,x_start,options_fminunc);
[x_2,fval_2] = fminsearch(f_con_xsi_mu,x_start,options_fminsearch);
if fval_1 < fval_2
  x_mle_constained_xsi_mu = x_1;
  fval = fval_1;
else
  x_mle_constained_xsi_mu = x_2;
  fval = fval_2;
end
llf_constrained_xsi_mu = -fval;
%  Compute adnjustment factor for test statistic
fname = [cv_dir 'cvx_test4_' str_cvs '_' num2str(T) '.csv'];
c_parm = readmatrix(fname);
adj = exp(c_parm(1) + c_parm(2)*xsi_mle + c_parm(3)*xsi_mle^2);
llf_dif = llf_unconstrained-llf_constrained_xsi_mu;
Stat_con_xsi_mu = adj*llf_dif; 

% Construct Confidence Interval for xsi;
% Form CI for xsi
% Step 1 Get Adjustment factor
% Adjustment factor
fname = [cv_dir 'cvx_test1_' str_cvs '_' num2str(T) '.csv'];
c_parm = readmatrix(fname);
adj = exp(c_parm(1) + c_parm(2)*xsi_mle + c_parm(3)*xsi_mle^2);
% Step 2: Input a grid of values for xsi
xsi_min = 0.5;
xsi_max = 1.3;
n_xsi = 500;
xsi_grid = linspace(xsi_min,xsi_max,n_xsi)';
% Step 3 construct log-likelihood value for each value in Xsi Grid
  llf_vec_con_xsi = NaN(n_xsi,1);
  llf_vec_con_xsi1 = NaN(n_xsi,1);
  x_start = [sigma_start mu_start]; % Starting Values
  for i = 1:n_xsi
    xsi_con = xsi_grid(i);
    f_con_xsi = @(x) minus_llf_exceedance_gev_con_xsi(x,Y_data,nobs_vec,tau_vec,xsi_con);
    [~,fval_1]=fminunc(f_con_xsi,x_start,options_fminunc);
    [~,fval_2]=fminsearch(f_con_xsi,x_start,options_fminsearch);
    if fval_1 < fval_2
      fval = fval_1;
    else
      fval = fval_2;
    end
    llf_vec_con_xsi(i)=-fval;
  end
  llf_dif = llf_unconstrained - llf_vec_con_xsi;
  stat = adj*llf_dif;
  ii = stat < 1.0;
  tmp = xsi_grid(ii==1);
  CI_xsi_lower = min(tmp);
  CI_xsi_upper = max(tmp);
  % Make sure these are not at bounds 
  if CI_xsi_lower == xsi_min
      CI_xsi_lower = NaN;
  end
  if CI_xsi_upper == xsi_max
      CI_xsi_upper = NaN;
  end

  % Construct a confidence interval for x_90 -- GEV 90th quantile
  % Get MLE
  x_90_mle = gevinv(0.9,xsi_mle,sigma_mle,mu_mle);
  % Adjustment factor
  fname = [cv_dir 'cvx_test2_' str_cvs '_' num2str(T) '.csv'];
  c_parm = readmatrix(fname);
  adj = exp(c_parm(1) + c_parm(2)*xsi_mle + c_parm(3)*xsi_mle^2);
  
  % Step 1: Input a grid of values 
  q_min = 4.0;
  q_max = 12.0;
  n_q = 500;
  q_grid = linspace(q_min,q_max,n_q)';
  % Step 2 construct log-likelihood value for each value in Xsi Grid
  llf_vec_con_q = NaN(n_q,1);
  % Starting value for MLE estimation
  x_start_mle = [xsi_mle sigma_mle];
  x_start = x_start_mle;

  for i = 1:n_q
    x_90_con = q_grid(i);
    f_con_x_90 = @(x) minus_llf_exceedance_gev_con_x_90(x,Y_data,nobs_vec,tau_vec,x_90_con);
    [x_1,fval_1]=fminunc(f_con_x_90,x_start,options_fminunc);
    [x_2,fval_2]=fminsearch(f_con_x_90,x_start,options_fminsearch);
    if fval_1 < fval_2
        fval = fval_1;
        x_1 = x_1;
    else
        fval = fval_2;
        x_1 = x_2;
    end
    llf_vec_con_q(i)=-fval;
    x_start = x_1;
 end

 llf_dif = llf_unconstrained - llf_vec_con_q;
  stat = adj*llf_dif;
  ii = stat < 1.0;
  tmp = q_grid(ii==1);
  CI_x90_lower = min(tmp);
  CI_x90_upper = max(tmp);
  % Make sure these are not at bounds 
  if CI_x90_lower == q_min
      CI_x90_lower = NaN;
  end
  if CI_x90_upper == q_max
      CI_x90_upper = NaN;
  end

% % Save Results to output file
outfname = [outdir application_str '_LikelihoodResults.txt'];
fid = fopen(outfname,'w');
fprintf(fid,'Data: %s\n',data_fname);
fprintf(fid,'T = %3i \n',size(Y_data,1));
fprintf(fid,'k = %3i \n',size(Y_data,2));
fprintf(fid,'MLEs with Hessian-based SEs \n');
fprintf(fid,'xsi_mle = %6.4f (%6.4f) \n',[xsi_mle se_xsi_mle]);
fprintf(fid,'sigma_mle = %6.4f (%6.4f) \n',[sigma_mle se_sigma_mle]);
fprintf(fid,'mu_mle = %6.4f (%6.4f) \n',[mu_mle se_mu_mle]);
fprintf(fid,'MLEs of lambda and Psi with tau = 1 \n');
[lambda_mle psi_mle] = GEV_to_PGP_Parms(xsi_mle,sigma_mle,mu_mle,1.0);
fprintf(fid,'      Lambda: %6.2f \n',lambda_mle);
fprintf(fid,'      Psi: %6.2f \n',psi_mle);
fprintf(fid,'MLEs 90th quantile of GEV: %6.2f \n',x_90_mle);
fprintf(fid,'95 percent CI for Xsi: %6.2f to %6.2f \n',[CI_xsi_lower CI_xsi_upper]);
fprintf(fid,'95 percent CI for 90th GEV quantile: %6.2f to %6.2f \n',[CI_x90_lower CI_x90_upper]);
fprintf(fid,'Some Test Statistics (5 percent CVs are 1.0) \n');
fprintf(fid,'   Test that mu = sigma/xsi: %6.2f \n',Stat_con_mu);
fprintf(fid,'   Test that mu = sigma/xsi and xsi = 1: %6.2f \n',Stat_con_xsi_mu);
fprintf(fid,'   Test for constant parameters (Nyblom)  %6.2f \n',adj_L_Stat);

fclose(fid);

