
% CODE FOR FIG S1

addpath('../')


%===========================================================================
% LOAD CLIMATE INDICES

% Load full ensembles for the nino and amo variables
efl='cesm_lme010_r12000_p998_state_2ozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:24:35.mat';
% LOAD JJA FOR AMO ANALYSIS
%efl='cesm_lme010_r12000_p998_state_2ozn_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:11:15.mat';

pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';
%load([pth,efl],'r_o','r_f','state_tp','xmeta','Xa_m','Xa_ens','reconYrs','p_yrs')
load([pth,efl])



% What years to look at?
analyrs1=1000:1925;
[~,~,ir1] = intersect(analyrs1,r_o:r_f); % indices of recon segment
% monthly indices
ovrlprng=rng_intr([r_o,r_f],[analyrs1(1),analyrs1(end)]);
mon_r=r_o:(1/12):(r_f+11/12);
irm1_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
irm1_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;



k=find(state_tp=='o');
xa_amo_ens=permute(squeeze(Xa_ens(xmeta{k}.id_X(1),:,:)),[2 1]);
% Ensemble mean of AMO
xa_amo_mn=mean(xa_amo_ens,2);
amo_off=mean(mean(xa_amo_ens(ir1,:),2));


k=find(state_tp=='2');
xa_gmt_ens=permute(squeeze(Xa_ens(xmeta{k}.id_X(1),:,:)),[2 1]);

k=find(state_tp=='n'); % specific nino variable
if strcmp(xmeta{k}.indinfo{3}.indxnm,'Nino 3.4'); i=3; else; error('Wrong Nino specification'); end
xa_nino_ens=reshape(permute(Xa_ens((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:,:),[1 3 2]),12*reconYrs,length(p_yrs));
xa_nino_mn=reshape(Xa_m((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:),12*reconYrs,1);
% Nino offset with respect to reference period
nino_off=mean(mean(xa_nino_ens(irm1_a:irm1_b,:),2));

% Annual nino from monthly data (leaves year ends not quite correct)
xa_nino_ann=mon2ann(xa_nino_mn,4,3);
% Annual nino from spatial field
%xa_nino_ann=nino(Xt2m,lat,lon,3,'n');

% Get annual Nino 3.4 value from annual spatial field
efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:33:07.mat';
load([pth,efl],'state_tp','xmeta','Xa_m')
k=find(state_tp=='t');
lat=xmeta{k}.lat;
lon=xmeta{k}.lon;
Xt2m=reshape(Xa_m(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:),length(lat),length(lon),reconYrs);

% Load mask for Chile/Argentina
%msk=load('ChlArg_mask.mat');
%msk=load('Pampas_mask.mat');
%lndmsk=msk.sftlf;
%lndmsk(lndmsk==0)=NaN;
% Land mask over Chile
msk=load('chile_mask.mat');
lndmsk=msk.chl;

% Latitude weighting
A=cosd(repmat(xmeta{k}.lat,[1 length(xmeta{k}.lon)]));



% Get the ensemble mean; identical to the original reconstruction
Xv_msk=bsxfun(@times,Xt2m,lndmsk);
t2m_casa=wmean_a(Xv_msk,A); % RECONSTRUCTION DROUGHT INDEX




% Drought atlases use JJA instead of annual
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:49:45.mat';
%pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';load([pth,efl])
%pth='/d2/nsteiger/output-da/hydroclimate/'; load([pth,efl])

% PDSI ONLY, ENSEMBLE RECONSTRUCTION
%efl='cesm_lme010_r12000_p998_state_q_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_06-Aug-2018_17:32:52.mat';
%pth='./output-recon/';load([pth,efl],'xmeta')


%=====================================
% LOAD ENSEMBLE MEGADROUGHT ESTIMATES

% Load both south american and north america drought indices
%load('svmd_ens_90pagree_sa.mat')
%load('svmd_ens_90pagree_pampas.mat')
%load('pampas_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('pampas_megaD_ens_noNAproxies.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('pampas_megaD_ens_noNASWproxies.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
load('chile_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')

xa_di_sa=Xa_di;
xa_di_ens_sa=Xa_di_ens;
ds_sa=svmd;
mdagr_sa=ens_drght;
sa_yrs=analyrs;

%load('svmd_ens_50pagree.mat')
%load('svmd_ens_90pagree.mat')
load('nasw_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
xa_di_na=Xa_di;
xa_di_ens_na=Xa_di_ens;
ds_na=svmd;
mdagr_na=ens_drght;
na_yrs=analyrs;


% Droughts in the real data mean
ds0_n=mdrought(xa_di_na(analyrs),1,'y');
ds0_s=mdrought(xa_di_sa(analyrs),1,'y');

% Droughts too frequent, throw out based on relative threshold
dcut=50;
np67=prctile(ds0_n(:,3),dcut);
rmv=find(ds0_n(:,3)>np67);
ds0_n(rmv,:)=[];

sp67=prctile(ds0_s(:,3),dcut);
rmv=find(ds0_s(:,3)>sp67);
ds0_s(rmv,:)=[];


%----------------------------
% Probabilites for real data
%----------------------------

% North America
dn0=ds0_n(:,1:2)';dn0=dn0(:);

% Fill out all years of droughts in vector
k=1;dn1=[];
while k<length(dn0)
   dn1=cat(2,dn1,dn0(k):dn0(k+1));
   k=k+2;
end

% South America
ds0=ds0_s(:,1:2)';ds0=ds0(:);

% Fill out all years of droughts in vector
k=1;ds1=[];
while k<length(ds0)
   ds1=cat(2,ds1,ds0(k):ds0(k+1));
   k=k+2;
end

% % 1/0 for kinds of drought
% dtm=zeros(length(analyrs),1); 
% dtm(dn1)=1;dtm(ds1)=1;ibth=intersect(dn1,ds1);dtm(ibth)=1;
% df0=find(dtm==0); % no drought will still be zero


% % Probabilities for differention possibilities
% pd_nthr0=length(df0)/length(analyrs);
% pd_na0=length(dn1)/length(analyrs);
% pd_sa0=length(ds1)/length(analyrs);
% pd_bh0=length(ibth)/length(analyrs);


% Figure out indices of neither, NA, SA, both
inthr=setdiff(1:length(analyrs),union(dn1,ds1)); % neither
[~,ina,isa]=setxor(dn1,ds1); % indices of NA only and SA only
ibth=intersect(dn1,ds1); % both

% Indices of alternatives
inotNA=setdiff(1:length(analyrs),dn1);
inotSA=setdiff(1:length(analyrs),ds1); 
inotNthr=setdiff(1:length(analyrs),inthr); 
inotBth=setdiff(1:length(analyrs),ibth); 




%===========================================================================
% LOAD FORCED CLIMATE RESPONSE ESTIMATES

% Full forcing
%load('../megad-forcing/output_twoboxCM_full_06-Aug-2018_15:09:13.mat','gmt_fr','myrs')
load('../ebm-held/output_twoboxCM_full_21-Sep-2018_13:29:49.mat','gmt_fr','myrs')
%load('../ebm-held/output_twoboxCM_full_SHvolcforc_eVolv2k_25-Nov-2019_15:50:31.mat','gmt_fr','myrs')
%load('../ebm-held/output_twoboxCM_full_eVolv2k_25-Nov-2019_16:02:15.mat','gmt_fr','myrs')
%load('../ebm-held/output_twoboxCM_full_eVolv2k_25-Nov-2019_16:02:15.mat','gmt_fr','myrs')
%load('../ebm-held/output_twoboxCM_full_SHvolcforcNotropics_26-Nov-2019_12:51:20.mat','gmt_fr','myrs')

% Plot mean forced response
gmt_est_mn=mean(gmt_fr,2);
[~,~,ig2] = intersect(analyrs1,myrs);% 1601-1925 reference 
%[~,~,ig3] = intersect(analyrs3,myrs); 
%[~,~,ig4] = intersect(analyrs4,myrs); 

% remove mean from ref
gmt_est_ens=gmt_fr-mean(gmt_est_mn(ig2));
gmt_est_f=gmt_est_mn-mean(gmt_est_mn(ig2));

% solar 
load('../ebm-held/output_twoboxCM_sol_21-Sep-2018_13:29:03.mat','gmt_fr')
%load('../megad-forcing/output_twoboxCM_solar_08-Aug-2018_13:47:28.mat','gmt_fr')
gmt_est_mn=mean(gmt_fr,2);
gmt_est_s=gmt_est_mn-mean(gmt_est_mn(ig2));

% volcanic 
%load('../megad-forcing/output_twoboxCM_volc_08-Aug-2018_13:38:26.mat','gmt_fr')
load('../ebm-held/output_twoboxCM_volc_21-Sep-2018_13:35:31.mat','gmt_fr')
gmt_est_mn=mean(gmt_fr,2);
gmt_est_v=gmt_est_mn-mean(gmt_est_mn(ig2));

% ghg
%load('../megad-forcing/output_twoboxCM_ghg_08-Aug-2018_14:05:19.mat','gmt_fr')
load('../ebm-held/output_twoboxCM_ghg_21-Sep-2018_13:32:35.mat','gmt_fr')
gmt_est_mn=mean(gmt_fr,2);
gmt_est_g=gmt_est_mn-mean(gmt_est_mn(ig2));

% ghg + solar
load('../ebm-held/output_twoboxCM_solghg_26-Sep-2018_14:06:47.mat','gmt_fr')
gmt_est_mn=mean(gmt_fr,2);
gmt_est_sg=gmt_est_mn-mean(gmt_est_mn(ig2));

% volc + solar
load('../ebm-held/output_twoboxCM_volcsol_26-Sep-2018_15:37:39.mat','gmt_fr')
gmt_est_mn=mean(gmt_fr,2);
gmt_est_vs=gmt_est_mn-mean(gmt_est_mn(ig2));

%===========================================================================


load distinct_colors_20.mat



%------------------------------------------
% LINEAR REGRESSION WITH AVERAGED DATA
%------------------------------------------


% Chile PDSI -----------------------------------------------------------

% Fit all varieties of models: What can explain PDSI in the American southwest?
%xm={1,2,3,4,5,[1,2],[1,3],[1,4],[1,5],[1,2,3],[1,2,4],[1,2,5]};
%c0={'Nino3.4','AMO','Forc','Vol','Sol','N+A','N+F','N+V','N+S','N+A+F','N+A+V','N+A+S'};
xm={1,2,[1,2]};


% Averaging length
%al=10;

bs=333;
%al=[8 9 10 11 12];
al=[9 10 11];
rsqrd=zeros(length(xm),bs,length(al));
aic_bs=zeros(length(xm),bs,length(al));
bic_bs=zeros(length(xm),bs,length(al));
dwp_avg4=zeros(length(xm),bs,length(al));
beta_bs=NaN(length(xm),4,bs,length(al));
p_bs=NaN(length(xm),4,bs,length(al));

k=1;
for j=al 
   % estimate NASW PDSI
   clear X4 y4
   %y4=detrend(ann2xt(Xa_di(ir1),j));
   y4=detrend(ann2xt(xa_di_sa(ir1),j));
   X4(:,1)=detrend(ann2xt(xa_nino_ann(ir1),j));
   %X4(:,2)=detrend(ann2xt(xa_amo_mn(ir1),j));
   X4(:,2)=detrend(ann2xt(gmt_est_f(ig2),j));
   %X4(:,2)=detrend(ann2xt(gmt_est_v(ig2),j));
   %X4(:,2)=detrend(ann2xt(gmt_est_s(ig2),j));
   %X4(:,2)=detrend(ann2xt(gmt_est_vs(ig2),j));

   %X4(:,3)=detrend(ann2xt(gmt_est_sg(ig4),j));
   %X4(:,3)=detrend(ann2xt(gmt_est_vs(ig4),j));
   %X4(:,4)=detrend(ann2xt(gmt_est_v(ig4),j)); 
   %X4(:,5)=detrend(ann2xt(gmt_est_s(ig4),j)); 
   % Standardize data to compare regression coefficients
   %X0 = bsxfun(@minus,X4,nanmean(X4,1));
   %X4 = bsxfun(@rdivide,X0,nanstd(X0,[],1));


   for i=1:length(xm)
      for ii=1:bs 
         % Bootstrap sample
         rsp=randsample(length(y4),length(y4),true);
         % Jacknife sample
%         rsp=randsample(length(y4),length(y4)-1);
         % ESTIMATE PDSI
	 mdl=fitlm(X4(rsp,xm{i}),y4(rsp));

	 rsqrd(i,ii,k)=mdl.Rsquared.Adjusted;
	 beta_bs(i,[1 (xm{i}+1)],ii,k)=mdl.Coefficients.Estimate;
	 p_bs(i,[1 (xm{i}+1)],ii,k)=mdl.Coefficients.pValue;
%	 dwp_avg4(i,ii,k)=dwtest(mdl);
	 aic_bs(i,ii,k)=mdl.ModelCriterion.AICc;   
	 bic_bs(i,ii,k)=mdl.ModelCriterion.BIC;   
       
      end
   end
   k=k+1;
end

r2_s=reshape(rsqrd,length(xm),bs*length(al));
% aic0=reshape(aic_bs,length(xm),bs*length(al));
% bic0=reshape(bic_bs,length(xm),bs*length(al));
% b0=reshape(beta_bs,length(xm),4,bs*length(al));
% p0=reshape(p_bs,length(xm),4,bs*length(al));
% 
% % Just take values for individual models
% clear b1 p1
% b1(:,1)=squeeze(b0(1,2,:));
% b1(:,2)=squeeze(b0(2,3,:));
% %b1(:,3)=squeeze(b0(3,4,:));
% %b1(:,4)=squeeze(b0(4,5,:));
% %b1(:,5)=squeeze(b0(5,6,:));
% 
% p1(:,1)=squeeze(p0(1,2,:));
% p1(:,2)=squeeze(p0(2,3,:));
% %p1(:,3)=squeeze(p0(3,4,:));
% %p1(:,4)=squeeze(p0(4,5,:));
% %p1(:,5)=squeeze(p0(5,6,:));



% NASW PDSI -----------------------------------------------------------

%xm={1,2,[1,2]};
xm={1,2,3,[1,2,3]}; % + AMO


% Averaging length
%al=10;

bs=333;
%al=[8 9 10 11 12];
al=[9 10 11];
rsqrd=zeros(length(xm),bs,length(al));
aic_bs=zeros(length(xm),bs,length(al));
bic_bs=zeros(length(xm),bs,length(al));
dwp_avg4=zeros(length(xm),bs,length(al));
beta_bs=NaN(length(xm),4,bs,length(al));
p_bs=NaN(length(xm),4,bs,length(al));

k=1;
for j=al 
   % estimate NASW PDSI
   clear X4 y4
   %y4=detrend(ann2xt(Xa_di(ir1),j));
   y4=detrend(ann2xt(xa_di_na(ir1),j));
   X4(:,1)=detrend(ann2xt(xa_nino_ann(ir1),j));
   X4(:,2)=detrend(ann2xt(xa_amo_mn(ir1),j));
   X4(:,3)=detrend(ann2xt(gmt_est_f(ig2),j));
   %X4(:,2)=detrend(ann2xt(gmt_est_v(ig2),j));
   %X4(:,2)=detrend(ann2xt(gmt_est_s(ig2),j));
   %X4(:,2)=detrend(ann2xt(gmt_est_vs(ig2),j));

   %X4(:,3)=detrend(ann2xt(gmt_est_sg(ig4),j));
   %X4(:,3)=detrend(ann2xt(gmt_est_vs(ig4),j));
   %X4(:,4)=detrend(ann2xt(gmt_est_v(ig4),j)); 
   %X4(:,5)=detrend(ann2xt(gmt_est_s(ig4),j)); 
   % Standardize data to compare regression coefficients
   %X0 = bsxfun(@minus,X4,nanmean(X4,1));
   %X4 = bsxfun(@rdivide,X0,nanstd(X0,[],1));


   for i=1:length(xm)
      for ii=1:bs 
         % Bootstrap sample
         rsp=randsample(length(y4),length(y4),true);
         % Jacknife sample
%         rsp=randsample(length(y4),length(y4)-1);
         % ESTIMATE PDSI
	 mdl=fitlm(X4(rsp,xm{i}),y4(rsp));

	 rsqrd(i,ii,k)=mdl.Rsquared.Adjusted;
	 beta_bs(i,[1 (xm{i}+1)],ii,k)=mdl.Coefficients.Estimate;
	 p_bs(i,[1 (xm{i}+1)],ii,k)=mdl.Coefficients.pValue;
%	 dwp_avg4(i,ii,k)=dwtest(mdl);
	 aic_bs(i,ii,k)=mdl.ModelCriterion.AICc;   
	 bic_bs(i,ii,k)=mdl.ModelCriterion.BIC;   
       
      end
   end
   k=k+1;
end

r2_n=reshape(rsqrd,length(xm),bs*length(al));

r2=cat(1,r2_n,r2_s);

%c0={'NINO3.4','Forc','N+F','NINO3.4','Forc','N+F'};
c0={'NINO3.4','AMO','Forc','N+A+F','NINO3.4','Forc','N+F'};
%c=categorical(c0,c0);

% Colors for boxplots
cm2=cat(1,clrs.green,clrs.blue,clrs.red,clrs.black,clrs.green,clrs.red,clrs.black);
%cm2=cat(1,clrs.green,clrs.red,clrs.black);
%cm2=cat(1,clrs.green,clrs.red,clrs.black,clrs.green,clrs.red,clrs.black);



% Distributions of the r2 values
figure
hold on
%hl=boxplot(r2',c0','symbol','','colors',cm2);
hl=boxplot(r2','Labels',c0,'symbol','','colors',cm2);
set(hl,'linew',2)
set(findobj(gcf,'-regexp','Tag','\w*Whisker'),'LineStyle','-')% whisker lines solid
ax=gca;
ylim([-0.025 0.65])
ax.YGrid = 'on';
%line('xdata',ax.XLim,'ydata',[0 0],'linestyle',':','color',clrs.grey,'linewidth',2)
ylabel('Adjusted r^2','fontsize',18)
set(gca,'FontSize',16,'XTickLabelRotation',45)
fig = gcf;fig.PaperUnits = 'inches';
%fig.PaperPosition = [0 2 4.5 4.5]; % [left bottom width height]
fig.PaperPosition = [0 2 5 4.5]; % add more vars

print(['/home/nsteiger/Dropbox/sa-phyda/sana_linear_r2_allModels_bootstrap.png'],'-dpng','-r300');


return

%------------------------------------------
% LOGISTIC REGRESSION ANALYSIS
%------------------------------------------

% get overlapping segements in time

% NOW FIGURE OUT WHICH YEARS ARE CONTINUOUS CHUNKS
drghts=struct([]); % initialize structure
negpos=ones(length(ir1),1);
negpos(ibth)=0; % variable is "both" years
% Find the timing of droughts
i=1;j=1;flag=0;
while i+1 <= length(ir1)
% Find the start of a drought
if negpos(i)==0
   d_start=i;
   i=i+1;
   % Keep going to the end of a drought
   while negpos(i)~=1
       i=i+1;
       if i+1 > length(ir1)
	   d_end=i-1;
	   drghts(j).strtend=[d_start d_end];
	   flag=1;
	   break
       end
   end
   if flag==1;break;end % get out of nested loops
   
   d_end=i-1;
   drghts(j).strtend=[d_start d_end];
   
   j=j+1;
end

i=i+1;
end
   
% Pull out years from structure
d=zeros(length(drghts),2);
for kk=1:length(drghts)
  d(kk,1:2)=drghts(kk).strtend;
end

return

% Bin the time series and do a logistric regression 
% Break up time series into single megadrought events and non-events

clear dst

dst=sort(d);

% Max allowable length
%ml=max(dst(:,2)-dst(:,1));
ml=ceil(mean(dst(:,2)-dst(:,1)));

% add fake end point for calculation to end
dst(end+1,1)=length(analyrs1)+1;
dst(end,2)=length(analyrs1)+11;

% 'dw' is similar in structure to 'ds'
dw=[];
y=1;j=1;
for i=1:length(dst)
   sl=length(y:dst(i,1));
   if sl<=ml
      dw(j,1)=y;
      dw(j,2)=dst(i,1)-1;
      j=j+1;
   elseif sl>ml
      % if it's too long, break it up into chuncks close to 'ml'
      dv0=ceil(sl/ml);
      if dv0==2;dv0=3;end
      dts0=round(linspace(y,dst(i,1),dv0));
      for jj=1:dv0-1
	 dw(j,1)=dts0(jj);
	 dw(j,2)=dts0(jj+1)-1;
         j=j+1;
      end
   end
   y=dst(i,2)+1;
end

% Select out segments for binning and regression
dwdd=sort(cat(1,dw,ds));
xD=zeros(length(dwdd),3);
w0=zeros(length(dwdd),1);
y0=zeros(length(dwdd),1);

% Binning 
for i=1:size(dwdd,1)
  [~,~,irab] = intersect((dwdd(i,1)+analyrs1(1)-1):(dwdd(i,2)+analyrs1(1)-1),r_o:r_f); 
  w0(i)=length(irab); % bin size 

% Mean of bins
 y0(i)=mean(Xa_di(irab));
 xD(i,1)=mean(xa_nino_ann(irab));
% xD(i,2)=mean(xa_amo_mn(irab));
 xD(i,2)=mean(gmt_est_f(irab)); % myrs and r_o start in the same year so this indexing is OK
 %xD(i,3)=mean(gmt_est_sg(irab)); % myrs and r_o start in the same year so this indexing is OK
 %xD(i,3)=mean(gmt_est_vs(irab)); % myrs and r_o start in the same year so this indexing is OK
 %xD(i,4)=mean(gmt_est_v(irab));
 %xD(i,5)=mean(gmt_est_s(irab));


end

%disp('Add gaussian random variable')
xD(:,3)=randn(size(dwdd,1),1);

% Test sensitivy of removing large volcanoes
%disp('Removing outlier volcanic segments')
xDnv=xD;
%volyr=[21,24,26,33];
%volyr=[24,33];
% Remove value that is more than three scaled median absolute deviations (MAD) away from the median
volyr=find(isoutlier(xD(:,3))==1); % i=[24,25,33,50], which are two MD and two non-MD 
%volyr=[];

xDnv(volyr,:)=NaN; % remove all volcanic year info
%xDnv(volyr,3)=NaN; % remove just forcing values
%xDnv(volyr,4)=NaN;


% Standardize data to compare regression coefficients
%Xd=zscore(xD);
X0 = bsxfun(@minus,xD,nanmean(xD,1));
Xd = bsxfun(@rdivide,X0,nanstd(X0,[],1));
X0 = bsxfun(@minus,xDnv,nanmean(xDnv,1));
Xdnv = bsxfun(@rdivide,X0,nanstd(X0,[],1));
% don't standardize (for getting actual coefficient values)
%Xd=xD;
%Xdnv=xDnv;

% Normalize weights?
wD=w0;
%wD=w0/max(w0);
%wD=w0/length(analyrs4)*100;
%wD=ones(size(w0));


% Which were the dry years?
yD=zeros(length(dwdd),1);
[~,id,~]=intersect(dwdd(:,1),ds(:,1));
yD(id)=1;
Yd=categorical(yD);

%figure;autocorr(yD);title('Is the drought time series autocorrelated?')


%--------------------
%    Regression
%--------------------

% Fit all varieties of models: What can explain PDSI in the American southwest?
% USE SAME AS ABOVE....
%xm={1,2,3,4,5,[1,2],[1,3],[1,4],[1,5],[1,2,3],[1,2,4],[1,2,5]};
%c0={'Nino3.4','AMO','Forc','Vol','Sol','N+A','N+F','N+V','N+S','N+A+F','N+A+V','N+A+S'};
%xm={1,2,3,[1,2],[1,3],[2,3],[1,2,3],4};
%c0={'Nino3.4','AMO','Forc','N+A','N+F','A+F','N+A+F','R'};
xm={1,2,[1,2],3};
c0={'Nino3.4','Forc','N+F','R'};
c=categorical(c0,c0);

% Jackknife sampling for regression coefficients
js=999;

aic_j=zeros(length(xm),js,1);
bic_j=zeros(length(xm),js,1);
b_j=NaN(length(xm),4,js);
%b_j=NaN(length(xm),5,js); % adding noise var
%p_j=NaN(length(xm),4,js);

aic_jnv=zeros(length(xm),js,1);
bic_jnv=zeros(length(xm),js,1);
b_jnv=NaN(length(xm),4,js);
%b_jnv=NaN(length(xm),5,js); % adding noise var
%p_jnv=NaN(length(xm),4,js);

% Weighted logistic regression
for i=1:js

   % Jacknife sample
   %rsp=randsample(length(Yd),length(Yd)-1);
   % Bootstrap
   rsp=randsample(length(Yd),length(Yd),true);

   % Random variable
   %Xd(:,4)=zscore(randn(size(dwdd,1),1));
   %Xdnv(:,4)=zscore(randn(size(dwdd,1),1));
   Xd(:,3)=zscore(randn(size(dwdd,1),1));
   Xdnv(:,3)=zscore(randn(size(dwdd,1),1));

   for ii=1:length(xm)
      % Regression
      mdl = fitglm(Xd(rsp,xm{ii}),Yd(rsp),'Distribution','binomial','Weights',wD(rsp));
      % Grab stats
      b_j(ii,[1 (xm{ii}+1)],i)=mdl.Coefficients.Estimate;
      %p_j(ii,[1 (xm{ii}+1)],i)=mdl.Coefficients.pValue;
      aic_j(ii,i)=mdl.ModelCriterion.AICc;
      bic_j(ii,i)=mdl.ModelCriterion.BIC;

      % Regression WITH VOLCANOES REMOVED
      mdl = fitglm(Xdnv(rsp,xm{ii}),Yd(rsp),'Distribution','binomial','Weights',wD(rsp));
      % Grab stats
      b_jnv(ii,[1 (xm{ii}+1)],i)=mdl.Coefficients.Estimate;
      %p_jnv(ii,[1 (xm{ii}+1)],i)=mdl.Coefficients.pValue;
      aic_jnv(ii,i)=mdl.ModelCriterion.AICc;
      bic_jnv(ii,i)=mdl.ModelCriterion.BIC;
   end
end


clear b1
b1(:,1)=squeeze(b_j(1,2,:));
b1(:,2)=squeeze(b_j(2,3,:));
%b1(:,3)=squeeze(b_j(3,4,:));
%b1(:,4)=squeeze(b_j(4,5,:));
%b1(:,5)=squeeze(b_j(5,6,:));

clear b01
b01(:,1)=squeeze(b_j(length(xm)-1,2,:));
b01(:,2)=squeeze(b_j(length(xm)-1,3,:));
%b01(:,3)=squeeze(b_j(7,4,:));

clear b2
b2(:,1)=squeeze(b_jnv(1,2,:));
b2(:,2)=squeeze(b_jnv(2,3,:));
%b2(:,3)=squeeze(b_jnv(3,4,:));

clear b02
b02(:,1)=squeeze(b_jnv(length(xm)-1,2,:));
b02(:,2)=squeeze(b_jnv(length(xm)-1,3,:));
%b02(:,3)=squeeze(b_jnv(7,4,:));

clear br
%br(:,1)=squeeze(b_j(length(xm),5,:));
%br(:,2)=squeeze(b_jnv(length(xm),5,:));
br(:,1)=squeeze(b_j(length(xm),4,:));
br(:,2)=squeeze(b_jnv(length(xm),4,:));


%clear p1
%p1(:,1)=squeeze(p_j(1,2,:));
%p1(:,2)=squeeze(p_j(2,3,:));
%p1(:,3)=squeeze(p_j(3,4,:));
%p1(:,4)=squeeze(p_j(4,5,:));
%p1(:,5)=squeeze(p_j(5,6,:));

% Model coefficients

% Restructure the data
bg(:,1)=b1(:,1);
bg(:,2)=b2(:,1);
bg(:,3)=b1(:,2);
bg(:,4)=b2(:,2);
%bg(:,5)=b1(:,3);
%bg(:,6)=b2(:,3);

% New categories
%c01={'Nino3.4-All','Nino3.4-NLV','AMO-All','AMO-NLV','Forc-All','Forc-NLV'};
c01={'Nino3.4-All','Nino3.4-NLV','Forc-All','Forc-NLV'};
c1=categorical(c01,c01);

% Color matrix
%cm=cat(1,clrs.green,clrs.green,clrs.blue,clrs.blue,clrs.red,clrs.red);
cm=cat(1,clrs.green,clrs.green,clrs.red,clrs.red);

figure
hold on
hl=boxplot(bg,c1,'Colors',cm,'Symbol','');
set(hl,'linew',2)
set(findobj(gcf,'-regexp','Tag','\w*Whisker'),'LineStyle','-')% whisker lines solid
ax=gca;
ax.YLim=[-4 2];
line('xdata',ax.XLim,'ydata',[0 0],'linestyle',':','color',clrs.grey,'linewidth',2)
ax.YGrid = 'on';
%ax.GridLineStyle = '--';
% Add percentage of values below:above 0
for j=1:size(bg,2);
   bpc=length(find(bg(:,j)<0))/length(bg(:,j))*100;
   text(j,ax.YLim(2)-0.5,[sprintf('%0.0f',bpc) ':' sprintf('%0.0f',100-bpc)],'HorizontalAlignment','center','fontsize',14)
end
ylabel('\beta_1','fontsize',16)
set(gca,'FontSize',14,'XTickLabelRotation',45)

fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 5 4]; % [left bottom width height]

%print(['figs/sa_logress_betas_bootstrap_V_NV.png'],'-dpng','-r300');


%figure
%boxplot(br)



% Best model via AIC and BIC?
[~,m_aic]=min(aic_j); % minimum of the AIC
[~,m_bic]=min(bic_j); % minimum of the AIC
hc_a=(histcounts(m_aic,0.5:(length(xm)+.5))./length(m_aic)).*100;
hc_b=(histcounts(m_bic,0.5:(length(xm)+.5))./length(m_bic)).*100;
hcab0=cat(1,hc_a,hc_b);

% NO VOLCANOES
[~,m_aic]=min(aic_jnv); % minimum of the AIC
[~,m_bic]=min(bic_jnv); % minimum of the AIC
hc_a=(histcounts(m_aic,0.5:(length(xm)+.5))./length(m_aic)).*100;
hc_b=(histcounts(m_bic,0.5:(length(xm)+.5))./length(m_bic)).*100;
hcab1=cat(1,hc_a,hc_b);

hcab=cat(1,hcab0,hcab1);

figure
b=bar(c,hcab','facecolor','flat','edgecolor','flat');
b(1).CData = clrs.white;
b(1).EdgeColor = clrs.grey;
b(1).LineWidth = 2;
b(2).CData = clrs.white;
b(2).EdgeColor = clrs.black;
b(2).LineWidth = 2;
b(3).CData = clrs.grey;
b(4).CData = clrs.black;
legend('AIC-All','BIC-All','AIC-NLV','BIC-NLV','location','northwest')
ylabel('Best model (%)','fontsize',16)
set(gca,'FontSize',14,'XTickLabelRotation',45)

fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 6 4]; % [left bottom width height]

%print(['figs/sa_logress_AICBIC_bootstrap_V_NV.png'],'-dpng','-r300');






































return

% % Model coefficients
% figure
% hold on
% boxplot(b1,c(1:3))
% ax=gca;
% %line('xdata',ax.XLim,'ydata',[0 0],'linestyle','--','color',clrs.grey,'linewidth',2)
% ax.YGrid = 'on';
% %ax.GridLineStyle = '--';
% ylabel('\beta_1','fontsize',16)
% set(gca,'FontSize',14,'XTickLabelRotation',45)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 3 4]; % [left bottom width height]
% 
% %print(['figs/linear_betas_jackknife.png'],'-dpng','-r300');
% %print(['figs/linear_betas_bootstrap_redux.png'],'-dpng','-r300');


% Model coefficients
% figure
% hold on
% boxplot(p1,c(1:5))
% ax=gca;
% line('xdata',ax.XLim,'ydata',[0.05 0.05],'linestyle','--','color',clrs.grey,'linewidth',2)
% ylabel('p-value','fontsize',16)
% set(gca,'FontSize',14,'XTickLabelRotation',45)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 4 4]; % [left bottom width height]

%print(['figs/linear_pval_jackknife.png'],'-dpng','-r300');
%print(['figs/linear_pval_bootstrap.png'],'-dpng','-r300');


% % Best model via AIC and BIC?
% [~,m_aic]=min(aic0); % minimum of the AIC
% [~,m_bic]=min(bic0); % minimum of the AIC
% hc_a=(histcounts(m_aic,0.5:(length(xm)+.5))./length(m_aic)).*100;
% hc_b=(histcounts(m_bic,0.5:(length(xm)+.5))./length(m_bic)).*100;
% hcab=cat(1,hc_a,hc_b);
% 
% figure
% bar(c,hcab')
% ylabel('Best model estimate (%)','fontsize',16)
% legend('AIC','BIC','location','southeast')
% set(gca,'FontSize',14,'XTickLabelRotation',45)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 5 4]; % [left bottom width height]
% 
% %print(['figs/linear_AICBIC_jackknife.png'],'-dpng','-r300');
% %print(['figs/linear_AICBIC_bootstrap_redux.png'],'-dpng','-r300');


%---------------------------------------------
% HOW WELL DO FORCINGS EXPLAIN NINO AND AMO?


xm={1,2,3};
c03={'Forc','Vol','Sol'};
c3=categorical(c03,c03);

% Averaging length
%al=10;

bs=333;
%al=[8 9 10 11 12];
al=[9 10 11];
rsqrd_n=zeros(length(xm),bs,length(al));
%rsqrd_a=zeros(length(xm),bs,length(al));
rsqrd_t=zeros(length(xm),bs,length(al));
beta_n=NaN(length(xm),4,bs,length(al));
%beta_a=NaN(length(xm),4,bs,length(al));
beta_t=NaN(length(xm),4,bs,length(al));

k=1;
for j=al 
   % estimate NASW PDSI
   clear X4 y4n y4a
   %y4n0=detrend(ann2xt(xa_nino_ann(ir1),j));
   %y4a0=detrend(ann2xt(xa_amo_mn(ir1),j));
   %y4t0=detrend(ann2xt(t2m_casa(ir1),j));
   y4n=detrend(ann2xt(xa_nino_ann(ir1),j));
  % y4a=detrend(ann2xt(xa_amo_mn(ir1),j));
   y4t=detrend(ann2xt(t2m_casa(ir1),j));

   X4(:,1)=detrend(ann2xt(gmt_est_f(ig4),j));
   %X4(:,1)=detrend(ann2xt(gmt_est_sg(ig4),j));
   %X4(:,1)=detrend(ann2xt(gmt_est_vs(ig4),j));

   X4(:,2)=detrend(ann2xt(gmt_est_v(ig4),j)); 
   X4(:,3)=detrend(ann2xt(gmt_est_s(ig4),j)); 
   % Standardize data to compare regression coefficients
   %y4n=(y4n0-nanmean(y4n0))./nanstd(y4n0);
   %y4a=(y4a0-nanmean(y4a0))./nanstd(y4a0);
   %y4t=(y4t0-nanmean(y4t0))./nanstd(y4t0);
   %X0 = bsxfun(@minus,X4,nanmean(X4,1));
   %X4 = bsxfun(@rdivide,X0,nanstd(X0,[],1));


   for i=1:length(xm)
      for ii=1:bs 
         % Bootstrap sample
         rsp=randsample(length(y4n),length(y4n),true);
         % Jacknife sample
%         rsp=randsample(length(y4),length(y4)-1);

         % ESTIMATE AMO
%	 mdl=fitlm(X4(rsp,xm{i}),y4a(rsp));
%	 rsqrd_a(i,ii,k)=mdl.Rsquared.Adjusted;
%	 beta_a(i,[1 (xm{i}+1)],ii,k)=mdl.Coefficients.Estimate;

         % ESTIMATE NINO
	 mdl=fitlm(X4(rsp,xm{i}),y4n(rsp));
	 rsqrd_n(i,ii,k)=mdl.Rsquared.Adjusted;
	 beta_n(i,[1 (xm{i}+1)],ii,k)=mdl.Coefficients.Estimate;

         % ESTIMATE LOCAL T
	 mdl=fitlm(X4(rsp,xm{i}),y4t(rsp));
	 rsqrd_t(i,ii,k)=mdl.Rsquared.Adjusted;
	 beta_t(i,[1 (xm{i}+1)],ii,k)=mdl.Coefficients.Estimate;

      end
   end
   k=k+1;
end

%r2a=reshape(rsqrd_a,length(xm),bs*length(al));
r2n=reshape(rsqrd_n,length(xm),bs*length(al));
r2t=reshape(rsqrd_t,length(xm),bs*length(al));
%b0a=reshape(beta_a,length(xm),4,bs*length(al));
b0n=reshape(beta_n,length(xm),4,bs*length(al));
b0t=reshape(beta_t,length(xm),4,bs*length(al));

% Just take values for individual models
clear b1a b1n b1t
%b1a(:,1)=squeeze(b0a(1,2,:));
%b1a(:,2)=squeeze(b0a(2,3,:));
%b1a(:,3)=squeeze(b0a(3,4,:));

b1n(:,1)=squeeze(b0n(1,2,:));
b1n(:,2)=squeeze(b0n(2,3,:));
b1n(:,3)=squeeze(b0n(3,4,:));

b1t(:,1)=squeeze(b0t(1,2,:));
b1t(:,2)=squeeze(b0t(2,3,:));
b1t(:,3)=squeeze(b0t(3,4,:));


% JUST LOOK AT FULL FORCING, NOT COMPONENTS
b4(:,1)=squeeze(b0n(1,2,:));
%b4(:,2)=squeeze(b0a(1,2,:));
b4(:,2)=squeeze(b0t(1,2,:));
rfrc(1,:)=r2n(1,:);
%rfrc(2,:)=r2a(1,:);
rfrc(2,:)=r2t(1,:);
%c02={'Nino3.4','AMO','T2m'};
c02={'Nino3.4','T2m'};
c2=categorical(c02,c02);

cm1=cat(1,clrs.green,clrs.black);

figure
hold on
hl=boxplot(rfrc',c2,'symbol','','colors',cm1);
set(hl,'linew',2)
set(findobj(gcf,'-regexp','Tag','\w*Whisker'),'LineStyle','-')% whisker lines solid
ax=gca;
line('xdata',ax.XLim,'ydata',[0 0],'linestyle',':','color',clrs.grey,'linewidth',2)
ylabel('Adjusted r^2','fontsize',16)
%ylim([-0.025 0.4])
set(gca,'FontSize',14,'XTickLabelRotation',45)
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 3 4]; % [left bottom width height]

%print(['figs/sa_linear_r2_forcNinoAmoT2m_bootstrap.png'],'-dpng','-r300');

% Model coefficients
figure
hl=boxplot(b4,c2,'symbol','','colors',cm1);
set(hl,'linew',2)
set(findobj(gcf,'-regexp','Tag','\w*Whisker'),'LineStyle','-')% whisker lines solid
ax=gca;
%ax.YLim=[-0.225 0.425];
line('xdata',ax.XLim,'ydata',[0 0],'linestyle',':','color',clrs.grey,'linewidth',2)
%ylim([-0. 0.4])
ax.YGrid = 'on';
% Add percentage of values below:above 0
for j=1:size(b4,2);
   bpc=length(find(b4(:,j)<0))/length(b4(:,j))*100;
   text(j,ax.YLim(2)-0.05,[sprintf('%0.0f',bpc) ':' sprintf('%0.0f',100-bpc)],'HorizontalAlignment','center','fontsize',14)
end

ylabel('\beta_1','fontsize',16)
set(gca,'FontSize',14,'XTickLabelRotation',45)

fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 3 4]; % [left bottom width height]

%print(['figs/sa_linear_beta_forcNinoAmoT2m_bootstrap.png'],'-dpng','-r300');


% % Distributions of the r2 values
% figure
% hold on
% boxplot(r2a',c)
% ax=gca;
% line('xdata',ax.XLim,'ydata',[0 0],'linestyle','--','color',clrs.grey,'linewidth',2)
% ylabel('Adjusted r^2','fontsize',16)
% ylim([-0.025 0.4])
% title('AMO')
% set(gca,'FontSize',14,'XTickLabelRotation',45)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 2.5 4]; % [left bottom width height]
% 
% %print(['figs/linear_forcamo_r2_jackknife.png'],'-dpng','-r300');
% %print(['figs/linear_forcamo_r2_bootstrap.png'],'-dpng','-r300');
% 
% % Distributions of the r2 values
% figure
% hold on
% boxplot(r2n',c)
% ax=gca;
% line('xdata',ax.XLim,'ydata',[0 0],'linestyle','--','color',clrs.grey,'linewidth',2)
% ylabel('Adjusted r^2','fontsize',16)
% ylim([-0.025 0.4])
% title('Nino3.4')
% set(gca,'FontSize',14,'XTickLabelRotation',45)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 2.5 4]; % [left bottom width height]
% 
% %print(['figs/linear_forcnino_r2_jackknife.png'],'-dpng','-r300');
% %print(['figs/linear_forcnino_r2_bootstrap.png'],'-dpng','-r300');
% 
% % Distributions of the r2 values
% figure
% hold on
% boxplot(r2t',c)
% ax=gca;
% line('xdata',ax.XLim,'ydata',[0 0],'linestyle','--','color',clrs.grey,'linewidth',2)
% ylabel('Adjusted r^2','fontsize',16)
% ylim([-0.025 0.4])
% title('T2m')
% set(gca,'FontSize',14,'XTickLabelRotation',45)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 2.5 4]; % [left bottom width height]
% 
% %print(['figs/linear_forct2m_r2_bootstrap.png'],'-dpng','-r300');
% 
% 




% 
% % Model coefficients
% figure
% hold on
% boxplot(b1n,c(1:3))
% ax=gca;
% %line('xdata',ax.XLim,'ydata',[0 0],'linestyle','--','color',clrs.grey,'linewidth',2)
% ax.YGrid = 'on';
% %ax.GridLineStyle = '--';
% ylabel('\beta_1','fontsize',16)
% ylim([-0.65 0.8])
% title('Nino3.4')
% set(gca,'FontSize',14,'XTickLabelRotation',45)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 2.5 4]; % [left bottom width height]
% 
% %print(['figs/linear_forcnino_betas_jackknife.png'],'-dpng','-r300');
% %print(['figs/linear_forcnino_betas_bootstrap.png'],'-dpng','-r300');
% 
% % Model coefficients
% figure
% hold on
% boxplot(b1t,c(1:3))
% ax=gca;
% %line('xdata',ax.XLim,'ydata',[0 0],'linestyle','--','color',clrs.grey,'linewidth',2)
% ax.YGrid = 'on';
% %ax.GridLineStyle = '--';
% ylabel('\beta_1','fontsize',16)
% ylim([-0.65 0.8])
% title('T2m')
% set(gca,'FontSize',14,'XTickLabelRotation',45)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 2.5 4]; % [left bottom width height]
% 
% %print(['figs/linear_forct2m_betas_bootstrap.png'],'-dpng','-r300');






%------------------------------------------
% LOGISTIC REGRESSION ANALYSIS
%------------------------------------------

% Bin the time series and do a logistric regression 
% Break up time series into single megadrought events and non-events

clear dst

dst=sort(ds);

% Max allowable length
%ml=max(dst(:,2)-dst(:,1));
ml=ceil(mean(dst(:,2)-dst(:,1)));

% add fake end point for calculation to end
dst(end+1,1)=length(analyrs4)+1;
dst(end,2)=length(analyrs4)+11;

% 'dw' is similar in structure to 'ds'
dw=[];
y=1;j=1;
for i=1:length(dst)
   sl=length(y:dst(i,1));
   if sl<=ml
      dw(j,1)=y;
      dw(j,2)=dst(i,1)-1;
      j=j+1;
   elseif sl>ml
      % if it's too long, break it up into chuncks close to 'ml'
      dv0=ceil(sl/ml);
      if dv0==2;dv0=3;end
      dts0=round(linspace(y,dst(i,1),dv0));
      for jj=1:dv0-1
	 dw(j,1)=dts0(jj);
	 dw(j,2)=dts0(jj+1)-1;
         j=j+1;
      end
   end
   y=dst(i,2)+1;
end

% Select out segments for binning and regression
dwdd=sort(cat(1,dw,ds));
xD=zeros(length(dwdd),3);
w0=zeros(length(dwdd),1);
y0=zeros(length(dwdd),1);

% Binning 
for i=1:size(dwdd,1)
  [~,~,irab] = intersect((dwdd(i,1)+analyrs1(1)-1):(dwdd(i,2)+analyrs1(1)-1),r_o:r_f); 
  w0(i)=length(irab); % bin size 

% Mean of bins
 y0(i)=mean(Xa_di(irab));
 xD(i,1)=mean(xa_nino_ann(irab));
% xD(i,2)=mean(xa_amo_mn(irab));
 xD(i,2)=mean(gmt_est_f(irab)); % myrs and r_o start in the same year so this indexing is OK
 %xD(i,3)=mean(gmt_est_sg(irab)); % myrs and r_o start in the same year so this indexing is OK
 %xD(i,3)=mean(gmt_est_vs(irab)); % myrs and r_o start in the same year so this indexing is OK
 %xD(i,4)=mean(gmt_est_v(irab));
 %xD(i,5)=mean(gmt_est_s(irab));


end

%disp('Add gaussian random variable')
xD(:,3)=randn(size(dwdd,1),1);

% Test sensitivy of removing large volcanoes
%disp('Removing outlier volcanic segments')
xDnv=xD;
%volyr=[21,24,26,33];
%volyr=[24,33];
% Remove value that is more than three scaled median absolute deviations (MAD) away from the median
volyr=find(isoutlier(xD(:,3))==1); % i=[24,25,33,50], which are two MD and two non-MD 
%volyr=[];

xDnv(volyr,:)=NaN; % remove all volcanic year info
%xDnv(volyr,3)=NaN; % remove just forcing values
%xDnv(volyr,4)=NaN;


% Standardize data to compare regression coefficients
%Xd=zscore(xD);
X0 = bsxfun(@minus,xD,nanmean(xD,1));
Xd = bsxfun(@rdivide,X0,nanstd(X0,[],1));
X0 = bsxfun(@minus,xDnv,nanmean(xDnv,1));
Xdnv = bsxfun(@rdivide,X0,nanstd(X0,[],1));
% don't standardize (for getting actual coefficient values)
%Xd=xD;
%Xdnv=xDnv;

% Normalize weights?
wD=w0;
%wD=w0/max(w0);
%wD=w0/length(analyrs4)*100;
%wD=ones(size(w0));


% Which were the dry years?
yD=zeros(length(dwdd),1);
[~,id,~]=intersect(dwdd(:,1),ds(:,1));
yD(id)=1;
Yd=categorical(yD);

%figure;autocorr(yD);title('Is the drought time series autocorrelated?')


%--------------------
%    Regression
%--------------------

% Fit all varieties of models: What can explain PDSI in the American southwest?
% USE SAME AS ABOVE....
%xm={1,2,3,4,5,[1,2],[1,3],[1,4],[1,5],[1,2,3],[1,2,4],[1,2,5]};
%c0={'Nino3.4','AMO','Forc','Vol','Sol','N+A','N+F','N+V','N+S','N+A+F','N+A+V','N+A+S'};
%xm={1,2,3,[1,2],[1,3],[2,3],[1,2,3],4};
%c0={'Nino3.4','AMO','Forc','N+A','N+F','A+F','N+A+F','R'};
xm={1,2,[1,2],3};
c0={'Nino3.4','Forc','N+F','R'};
c=categorical(c0,c0);

% Jackknife sampling for regression coefficients
js=999;

aic_j=zeros(length(xm),js,1);
bic_j=zeros(length(xm),js,1);
b_j=NaN(length(xm),4,js);
%b_j=NaN(length(xm),5,js); % adding noise var
%p_j=NaN(length(xm),4,js);

aic_jnv=zeros(length(xm),js,1);
bic_jnv=zeros(length(xm),js,1);
b_jnv=NaN(length(xm),4,js);
%b_jnv=NaN(length(xm),5,js); % adding noise var
%p_jnv=NaN(length(xm),4,js);

% Weighted logistic regression
for i=1:js

   % Jacknife sample
   %rsp=randsample(length(Yd),length(Yd)-1);
   % Bootstrap
   rsp=randsample(length(Yd),length(Yd),true);

   % Random variable
   %Xd(:,4)=zscore(randn(size(dwdd,1),1));
   %Xdnv(:,4)=zscore(randn(size(dwdd,1),1));
   Xd(:,3)=zscore(randn(size(dwdd,1),1));
   Xdnv(:,3)=zscore(randn(size(dwdd,1),1));

   for ii=1:length(xm)
      % Regression
      mdl = fitglm(Xd(rsp,xm{ii}),Yd(rsp),'Distribution','binomial','Weights',wD(rsp));
      % Grab stats
      b_j(ii,[1 (xm{ii}+1)],i)=mdl.Coefficients.Estimate;
      %p_j(ii,[1 (xm{ii}+1)],i)=mdl.Coefficients.pValue;
      aic_j(ii,i)=mdl.ModelCriterion.AICc;
      bic_j(ii,i)=mdl.ModelCriterion.BIC;

      % Regression WITH VOLCANOES REMOVED
      mdl = fitglm(Xdnv(rsp,xm{ii}),Yd(rsp),'Distribution','binomial','Weights',wD(rsp));
      % Grab stats
      b_jnv(ii,[1 (xm{ii}+1)],i)=mdl.Coefficients.Estimate;
      %p_jnv(ii,[1 (xm{ii}+1)],i)=mdl.Coefficients.pValue;
      aic_jnv(ii,i)=mdl.ModelCriterion.AICc;
      bic_jnv(ii,i)=mdl.ModelCriterion.BIC;
   end
end


clear b1
b1(:,1)=squeeze(b_j(1,2,:));
b1(:,2)=squeeze(b_j(2,3,:));
%b1(:,3)=squeeze(b_j(3,4,:));
%b1(:,4)=squeeze(b_j(4,5,:));
%b1(:,5)=squeeze(b_j(5,6,:));

clear b01
b01(:,1)=squeeze(b_j(length(xm)-1,2,:));
b01(:,2)=squeeze(b_j(length(xm)-1,3,:));
%b01(:,3)=squeeze(b_j(7,4,:));

clear b2
b2(:,1)=squeeze(b_jnv(1,2,:));
b2(:,2)=squeeze(b_jnv(2,3,:));
%b2(:,3)=squeeze(b_jnv(3,4,:));

clear b02
b02(:,1)=squeeze(b_jnv(length(xm)-1,2,:));
b02(:,2)=squeeze(b_jnv(length(xm)-1,3,:));
%b02(:,3)=squeeze(b_jnv(7,4,:));

clear br
%br(:,1)=squeeze(b_j(length(xm),5,:));
%br(:,2)=squeeze(b_jnv(length(xm),5,:));
br(:,1)=squeeze(b_j(length(xm),4,:));
br(:,2)=squeeze(b_jnv(length(xm),4,:));


%clear p1
%p1(:,1)=squeeze(p_j(1,2,:));
%p1(:,2)=squeeze(p_j(2,3,:));
%p1(:,3)=squeeze(p_j(3,4,:));
%p1(:,4)=squeeze(p_j(4,5,:));
%p1(:,5)=squeeze(p_j(5,6,:));

% Model coefficients

% Restructure the data
bg(:,1)=b1(:,1);
bg(:,2)=b2(:,1);
bg(:,3)=b1(:,2);
bg(:,4)=b2(:,2);
%bg(:,5)=b1(:,3);
%bg(:,6)=b2(:,3);

% New categories
%c01={'Nino3.4-All','Nino3.4-NLV','AMO-All','AMO-NLV','Forc-All','Forc-NLV'};
c01={'Nino3.4-All','Nino3.4-NLV','Forc-All','Forc-NLV'};
c1=categorical(c01,c01);

% Color matrix
%cm=cat(1,clrs.green,clrs.green,clrs.blue,clrs.blue,clrs.red,clrs.red);
cm=cat(1,clrs.green,clrs.green,clrs.red,clrs.red);

figure
hold on
hl=boxplot(bg,c1,'Colors',cm,'Symbol','');
set(hl,'linew',2)
set(findobj(gcf,'-regexp','Tag','\w*Whisker'),'LineStyle','-')% whisker lines solid
ax=gca;
ax.YLim=[-4 2];
line('xdata',ax.XLim,'ydata',[0 0],'linestyle',':','color',clrs.grey,'linewidth',2)
ax.YGrid = 'on';
%ax.GridLineStyle = '--';
% Add percentage of values below:above 0
for j=1:size(bg,2);
   bpc=length(find(bg(:,j)<0))/length(bg(:,j))*100;
   text(j,ax.YLim(2)-0.5,[sprintf('%0.0f',bpc) ':' sprintf('%0.0f',100-bpc)],'HorizontalAlignment','center','fontsize',14)
end
ylabel('\beta_1','fontsize',16)
set(gca,'FontSize',14,'XTickLabelRotation',45)

fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 5 4]; % [left bottom width height]

%print(['figs/sa_logress_betas_bootstrap_V_NV.png'],'-dpng','-r300');


%figure
%boxplot(br)



% Best model via AIC and BIC?
[~,m_aic]=min(aic_j); % minimum of the AIC
[~,m_bic]=min(bic_j); % minimum of the AIC
hc_a=(histcounts(m_aic,0.5:(length(xm)+.5))./length(m_aic)).*100;
hc_b=(histcounts(m_bic,0.5:(length(xm)+.5))./length(m_bic)).*100;
hcab0=cat(1,hc_a,hc_b);

% NO VOLCANOES
[~,m_aic]=min(aic_jnv); % minimum of the AIC
[~,m_bic]=min(bic_jnv); % minimum of the AIC
hc_a=(histcounts(m_aic,0.5:(length(xm)+.5))./length(m_aic)).*100;
hc_b=(histcounts(m_bic,0.5:(length(xm)+.5))./length(m_bic)).*100;
hcab1=cat(1,hc_a,hc_b);

hcab=cat(1,hcab0,hcab1);

figure
b=bar(c,hcab','facecolor','flat','edgecolor','flat');
b(1).CData = clrs.white;
b(1).EdgeColor = clrs.grey;
b(1).LineWidth = 2;
b(2).CData = clrs.white;
b(2).EdgeColor = clrs.black;
b(2).LineWidth = 2;
b(3).CData = clrs.grey;
b(4).CData = clrs.black;
legend('AIC-All','BIC-All','AIC-NLV','BIC-NLV','location','northwest')
ylabel('Best model (%)','fontsize',16)
set(gca,'FontSize',14,'XTickLabelRotation',45)

fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 6 4]; % [left bottom width height]

%print(['figs/sa_logress_AICBIC_bootstrap_V_NV.png'],'-dpng','-r300');














% %---------------------------------------------
% % HOW WELL DO FORCINGS EXPLAIN NINO AND AMO?
% 
% 
% %xm={1,2,3};
% %c0={'Forc','Vol','Sol'};
% c0={'Forc'};
% c=categorical(c0,c0);
% 
% % Averaging length
% %al=10;
% 
% bs=200;
% %al=[8 9 10 11 12];
% al=[9 10 11];
% rsqrd_n=zeros(bs,length(al));
% rsqrd_a=zeros(bs,length(al));
% beta_n=NaN(2,bs,length(al));
% beta_a=NaN(2,bs,length(al));
% 
% k=1;
% for j=al 
%    % estimate NASW PDSI
%    clear X4 y4n y4a
%    y4n0=detrend(ann2xt(xa_nino_ann(ir4),j));
%    y4a0=detrend(ann2xt(xa_amo_mn(ir4),j));
%    X4=detrend(ann2xt(gmt_est_f(ig4),j));
%    %X4(:,2)=detrend(ann2xt(gmt_est_v(ig4),j)); 
%    %X4(:,3)=detrend(ann2xt(gmt_est_s(ig4),j)); 
%    % Standardize data to compare regression coefficients
%    y4n=(y4n0-nanmean(y4n0))./nanstd(y4n0);
%    y4a=(y4a0-nanmean(y4a0))./nanstd(y4a0);
%    X0 = bsxfun(@minus,X4,nanmean(X4,1));
%    X4 = bsxfun(@rdivide,X0,nanstd(X0,[],1));
% 
% 
%    for ii=1:bs 
%       % Bootstrap sample
%       rsp=randsample(length(y4n),length(y4n),true);
%       % Jacknife sample
% %         rsp=randsample(length(y4),length(y4)-1);
% 
%       % ESTIMATE AMO
%       mdl=fitlm(X4(rsp),y4a(rsp));
%       rsqrd_a(ii,k)=mdl.Rsquared.Adjusted;
%       beta_a(:,ii,k)=mdl.Coefficients.Estimate;
% 
%       % ESTIMATE NINO
%       mdl=fitlm(X4(rsp,xm{i}),y4n(rsp));
%       rsqrd_n(ii,k)=mdl.Rsquared.Adjusted;
%       beta_n(:,ii,k)=mdl.Coefficients.Estimate;
% 
%    end
%    k=k+1;
% end
% 
% r2a=reshape(rsqrd_a,bs*length(al));
% r2n=reshape(rsqrd_n,bs*length(al));
% b0a=reshape(beta_a,4,bs*length(al));
% b0n=reshape(beta_n,4,bs*length(al));
% 
% % Just take values for individual models
% clear b1a b1n
% b1a(:,1)=squeeze(b0a(1,2,:));
% %b1a(:,2)=squeeze(b0a(2,3,:));
% %b1a(:,3)=squeeze(b0a(3,4,:));
% 
% b1n(:,1)=squeeze(b0n(1,2,:));
% %b1n(:,2)=squeeze(b0n(2,3,:));
% %b1n(:,3)=squeeze(b0n(3,4,:));
% 
% 
% 
% % Distributions of the r2 values
% figure
% hold on
% boxplot(r2a',c)
% plot(1:length(c),zeros(length(c),1),'--','color',clrs.grey,'linewidth',2)
% ylabel('Adjusted r^2','fontsize',16)
% title('AMO')
% set(gca,'FontSize',14,'XTickLabelRotation',45)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 5 4]; % [left bottom width height]
% 
% %print(['figs/linear_forcamo_r2_jackknife.png'],'-dpng','-r300');
% print(['figs/linear_forcamo_r2_bootstrap.png'],'-dpng','-r300');
% 
% % Distributions of the r2 values
% figure
% hold on
% boxplot(r2n',c)
% plot(1:length(c),zeros(length(c),1),'--','color',clrs.grey,'linewidth',2)
% ylabel('Adjusted r^2','fontsize',16)
% title('Nino3.4')
% set(gca,'FontSize',14,'XTickLabelRotation',45)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 5 4]; % [left bottom width height]
% 
% %print(['figs/linear_forcnino_r2_jackknife.png'],'-dpng','-r300');
% print(['figs/linear_forcnino_r2_bootstrap.png'],'-dpng','-r300');
% 
% % Model coefficients
% figure
% hold on
% boxplot(b1a,c(1:3))
% ax=gca;
% %line('xdata',ax.XLim,'ydata',[0 0],'linestyle','--','color',clrs.grey,'linewidth',2)
% ax.YGrid = 'on';
% %ax.GridLineStyle = '--';
% ylabel('\beta_1','fontsize',16)
% title('AMO')
% set(gca,'FontSize',14,'XTickLabelRotation',45)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 4 4]; % [left bottom width height]
% 
% %print(['figs/linear_forcamo_betas_jackknife.png'],'-dpng','-r300');
% print(['figs/linear_forcamo_betas_bootstrap.png'],'-dpng','-r300');
% 
% % Model coefficients
% figure
% hold on
% boxplot(b1n,c(1:3))
% ax=gca;
% %line('xdata',ax.XLim,'ydata',[0 0],'linestyle','--','color',clrs.grey,'linewidth',2)
% ax.YGrid = 'on';
% %ax.GridLineStyle = '--';
% ylabel('\beta_1','fontsize',16)
% title('Nino3.4')
% set(gca,'FontSize',14,'XTickLabelRotation',45)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 4 4]; % [left bottom width height]
% 
% %print(['figs/linear_forcnino_betas_jackknife.png'],'-dpng','-r300');
% print(['figs/linear_forcnino_betas_bootstrap.png'],'-dpng','-r300');


%---------------------------------------
% DIFFERENT TIME PERIODS FOR REGRESSION

% % estimate NASW PDSI
% y1=detrend(ann2xt(Xa_di(ir1),al));
% X1(:,1)=detrend(ann2xt(xa_nino_ann(ir1),al));
% X1(:,2)=detrend(ann2xt(xa_amo_mn(ir1),al));
% X1(:,3)=detrend(ann2xt(gmt_est_f(ig1),al)); % full forcing
% X1(:,4)=detrend(ann2xt(gmt_est_v(ig1),al)); % full forcing
% X1(:,5)=detrend(ann2xt(gmt_est_s(ig1),al)); % full forcing
% 
% % estimate NASW PDSI
% y2=detrend(ann2xt(Xa_di(ir2),al));
% X2(:,1)=detrend(ann2xt(xa_nino_ann(ir2),al));
% X2(:,2)=detrend(ann2xt(xa_amo_mn(ir2),al));
% X2(:,3)=detrend(ann2xt(gmt_est_f(ig2),al)); % full forcing
% X2(:,4)=detrend(ann2xt(gmt_est_v(ig2),al)); % full forcing
% X2(:,5)=detrend(ann2xt(gmt_est_s(ig2),al)); % full forcing
% 
% % estimate NASW PDSI
% y4=detrend(ann2xt(Xa_di(ir4),al));
% X4(:,1)=detrend(ann2xt(xa_nino_ann(ir4),al));
% X4(:,2)=detrend(ann2xt(xa_amo_mn(ir4),al));
% X4(:,3)=detrend(ann2xt(gmt_est_f(ig4),al)); % full forcing
% X4(:,4)=detrend(ann2xt(gmt_est_v(ig4),al)); % full forcing
% X4(:,5)=detrend(ann2xt(gmt_est_s(ig4),al)); % full forcing
% 
% % y4o=ann2xt(Xa_di(ir4),al);
% % X4o(:,1)=ann2xt(xa_nino_ann(ir4),al);
% % X4o(:,2)=ann2xt(xa_amo_mn(ir4),al);
% % X4o(:,3)=ann2xt(gmt_est_f(ig4),al); % full forcing
% % X4o(:,4)=ann2xt(gmt_est_v(ig4),al); % full forcing
% % X4o(:,5)=ann2xt(gmt_est_s(ig4),al); % full forcing
% 
% rsqrd_avg1=zeros(length(xm),1);
% rsqrd_avg2=zeros(length(xm),1);
% rsqrd=zeros(length(xm),1);
% %rsqrdo=zeros(length(xm),1);
% aic_bs_avg1=zeros(length(xm),1);
% aic_bs_avg2=zeros(length(xm),1);
% aic_bs=zeros(length(xm),1);
% dwp_avg1=zeros(length(xm),1);
% dwp_avg2=zeros(length(xm),1);
% dwp_avg4=zeros(length(xm),1);
% %dwp_avg4o=zeros(length(xm),1);
% for i=1:length(xm);
%    mdl=fitlm(X1(:,xm{i}),y1);
%    rsqrd_avg1(i)=mdl.Rsquared.Adjusted;
%    dwp_avg1(i)=dwtest(mdl);
%    aic_bs_avg1(i)=mdl.ModelCriterion.AICc;   
% 
%    mdl=fitlm(X2(:,xm{i}),y2);
%    rsqrd_avg2(i)=mdl.Rsquared.Adjusted;
%    dwp_avg2(i)=dwtest(mdl);
%    aic_bs_avg2(i)=mdl.ModelCriterion.AICc;   
% 
%    mdl=fitlm(X4(:,xm{i}),y4);
%    rsqrd(i)=mdl.Rsquared.Adjusted;
%    dwp_avg4(i)=dwtest(mdl);
%    aic_bs(i)=mdl.ModelCriterion.AICc;   
% 
%    %mdl=fitlm(X4o(:,xm{i}),y4o);
%    %rsqrdo(i)=mdl.Rsquared.Adjusted;
%    %dwp_avg4o(i)=dwtest(mdl);
% 
% end
% 
% % horizontal bar plot
% figure
% % ax1=subplot(3,1,1);
% % barh(c,rsqrd_avg1)
% % title('800-1600')
% % ax1.XGrid='on';
% % xlabel(['Adjusted r^2'])
% % 
% % ax2=subplot(3,1,2);
% % barh(c,rsqrd_avg2)
% % title('1601-1925')
% % ax2.XGrid='on';
% % xlabel(['Adjusted r^2'])
% 
% %ax3=subplot(3,1,3);
% barh(c,rsqrd)
% title('800-1925')
% %ax3.XGrid='on';
% xlabel(['Adjusted r^2'])
% 
% %linkaxes([ax3,ax2,ax1],'x');
% 
% %set(gca,'fontsize',16)
% 
% %fig = gcf;fig.PaperUnits = 'inches';
% %fig.PaperPosition = [0 2 5 7.5]; % [left bottom width height]
% 
% %print(['figs/linearfit_comparison_r2_10yr.png'],'-dpng','-r300');
% 




%------------------------------------
% LINEAR REGRESSION WITH ANNUAL DATA
%------------------------------------


% % estimate NASW PDSI
% y1=detrend(Xa_di(ir1));
% X1(:,1)=detrend(xa_nino_ann(ir1));
% X1(:,2)=detrend(xa_amo_mn(ir1));
% X1(:,3)=detrend(gmt_est_f(ig1)); % full forcing
% X1(:,4)=detrend(gmt_est_v(ig1)); % full forcing
% X1(:,5)=detrend(gmt_est_s(ig1)); % full forcing
% 
% % estimate NASW PDSI
% y2=detrend(Xa_di(ir2));
% X2(:,1)=detrend(xa_nino_ann(ir2));
% X2(:,2)=detrend(xa_amo_mn(ir2));
% X2(:,3)=detrend(gmt_est_f(ig2)); % full forcing
% X2(:,4)=detrend(gmt_est_v(ig2)); % full forcing
% X2(:,5)=detrend(gmt_est_s(ig2)); % full forcing
% 
% % estimate NASW PDSI
% y4=detrend(Xa_di(ir4));
% X4(:,1)=detrend(xa_nino_ann(ir4));
% X4(:,2)=detrend(xa_amo_mn(ir4));
% X4(:,3)=detrend(gmt_est_f(ig4)); % full forcing
% X4(:,4)=detrend(gmt_est_v(ig4)); % full forcing
% X4(:,5)=detrend(gmt_est_s(ig4)); % full forcing
% 
% % ACCESS LINEAR MODEL FIT OUTPUT:
% % http://www.mathworks.com/help/stats/linearmodel-class.html
% 
% rsqrd_ann1=zeros(length(xm),1);
% rsqrd_ann2=zeros(length(xm),1);
% rsqrd_ann4=zeros(length(xm),1);
% aic_bs_ann1=zeros(length(xm),1);
% aic_bs_ann2=zeros(length(xm),1);
% aic_bs_ann4=zeros(length(xm),1);
% dwp_ann1=zeros(length(xm),1);
% dwp_ann2=zeros(length(xm),1);
% dwp_ann4=zeros(length(xm),1);
% for i=1:length(xm);
%    mdl=fitlm(X1(:,xm{i}),y1);
%    rsqrd_ann1(i)=mdl.Rsquared.Adjusted;
%    aic_bs_ann1(i)=mdl.ModelCriterion.AICc;   
%    dwp_ann1(i)=dwtest(mdl);
% 
%    mdl=fitlm(X2(:,xm{i}),y2);
%    rsqrd_ann2(i)=mdl.Rsquared.Adjusted;
%    aic_bs_ann2(i)=mdl.ModelCriterion.AICc;   
%    dwp_ann2(i)=dwtest(mdl);
%    
%    mdl=fitlm(X4(:,xm{i}),y4);
%    rsqrd_ann4(i)=mdl.Rsquared.Adjusted;
%    aic_bs_ann4(i)=mdl.ModelCriterion.AICc;   
%    dwp_ann4(i)=dwtest(mdl);
% end
% 
% % horizontal bar plot
% figure
% ax1=subplot(3,1,1);
% barh(c,rsqrd_ann1)
% title('800-1600')
% ax1.XGrid='on';
% xlabel(['Adjusted r^2'])
% 
% ax2=subplot(3,1,2);
% barh(c,rsqrd_ann2)
% title('1601-1925')
% ax2.XGrid='on';
% xlabel(['Adjusted r^2'])
% 
% ax3=subplot(3,1,3);
% barh(c,rsqrd_ann4)
% title('800-1925')
% ax3.XGrid='on';
% xlabel(['Adjusted r^2'])
% 
% linkaxes([ax3,ax2,ax1],'x');
% 
% %set(gca,'fontsize',16)
% 
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 5 7.5]; % [left bottom width height]
% 
% %print(['figs/linearfit_comparison_r2_ann.png'],'-dpng','-r300');
% 
% 
% % horizontal bar plot
% figure
% ax1=subplot(3,1,1);
% barh(c,1e-3.*aic_bs_ann1)
% title('800-1600')
% ax1.XGrid='on';
% xlabel(['AIC'])
% 
% ax2=subplot(3,1,2);
% barh(c,1e-3.*aic_bs_ann2)
% title('1601-1925')
% ax2.XGrid='on';
% xlabel(['AIC'])
% 
% ax3=subplot(3,1,3);
% barh(c,1e-3.*aic_bs_ann4)
% title('800-1925')
% ax3.XGrid='on';
% xlabel(['AIC'])
% 
% linkaxes([ax3,ax2,ax1],'x');
% 
% %set(gca,'fontsize',16)
% 
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 5 7.5]; % [left bottom width height]
% 
% %print(['figs/linearfit_comparison_r2_ann.png'],'-dpng','-r300');


% is proportional to the probability that the ith model minimizes the (estimated) information loss
% exp((AICmin − AICi)/2)


% % Can any forcings explain Nino or AMO?
% xm={3,4,5};
% rsqrd_n=zeros(length(xm),1);
% rsqrd_a=zeros(length(xm),1);
% for i=1:length(xm);
%    mdl=fitlm(X(:,xm{i}),y2);rsqrd_n(i)=mdl.Rsquared.Adjusted;
%    mdl=fitlm(X(:,xm{i}),y3);rsqrd_a(i)=mdl.Rsquared.Adjusted;
% end
% 
% 
% % horizontal bar plot
% c0={'Forc','Vol','Sol'};
% c=categorical(c0,c0);
% figure
% barh(c,rsqrd_n)
% title('Forcing explain Nino3.4?')
% xlabel(['Adjusted r^2'],'fontsize',18)
% set(gca,'fontsize',16)
% 
% %print(['figs/linearfit_nino_forcings_r2.png'],'-dpng','-r300');
% 
% figure
% barh(c,rsqrd_a)
% title('Forcing explain AMO?')
% xlabel(['Adjusted r^2'],'fontsize',18)
% set(gca,'fontsize',16)
% 
% %print(['figs/linearfit_amo_forcings_r2.png'],'-dpng','-r300');










% % Averaged PDSI
% y0=Xa_di(ir4)-mean(Xa_di(ir2));
% 
% % Break up time series by 'bl' year segments; find how many are drought years and bin accordingly
% bl=10; % bin length; 21 = smallest bin length with no autocorrelation in the period 800-1600
% segs=fix(length(y0)/bl);
% fn=length(y0)-mod(length(y0),bl);
% 
% ynD=sum(reshape(y0(1:fn),bl,segs))./bl;
% 
% % Mean climate and forcing values in each bin
% X=zeros(segs,6);
% X0=xa_nino_ann(ir4)-mean(xa_nino_ann(ir2));
% X(1:segs,1)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,1)=min(reshape(X0(1:fn),bl,segs));
% 
% X0=xa_amo_mn(ir4)-mean(xa_amo_mn(ir2));
% X(1:segs,2)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,2)=sum(reshape(X0(1:fn),bl,segs));
% 
% X0=gmt_est_f(ig4)-mean(gmt_est_f(ig2));
% X(1:segs,3)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,3)=min(reshape(X0(1:fn),bl,segs));
% 
% X0=gmt_est_v(ig4)-mean(gmt_est_v(ig2));
% X(1:segs,4)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,4)=sum(reshape(X0(1:fn),bl,segs));
% 
% X0=gmt_est_s(ig4)-mean(gmt_est_s(ig2));
% X(1:segs,5)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,5)=sum(reshape(X0(1:fn),bl,segs));
% 
% X0=gmt_est_g(ig4)-mean(gmt_est_g(ig2));
% X(1:segs,6)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,6)=sum(reshape(X0(1:fn),bl,segs));
% 
% 
% 
% % Fit all varieties of models: What can explain PDSI in the American southwest?
% xm={1,2,3,4,5,6,[1,2],[1,3],[1,4],[1,5],[1,6],[1,2,3],[1,2,4]};
% rsqrd=zeros(length(xm),1);
% for i=1:length(xm);
%    mdl=fitlm(X(:,xm{i}),ynD);rsqrd(i)=mdl.Rsquared.Adjusted;
% end
% 
% % horizontal bar plot
% c0={'Nino','AMO','Forc','Vol','Sol','GHG','N+A','N+F','N+V','N+S','N+G','N+A+F','N+A+V'};
% c=categorical(c0,c0);
% figure
% barh(c,rsqrd)
% grid on
% xlabel(['Adjusted r^2'],'fontsize',16)
% set(gca,'fontsize',14)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]
% title(['Linear regression, ' num2str(bl) ' yr binned data'])
% 
% 
% 
% % Can any forcings explain Nino or AMO?
% xm={3,4,5,6};
% rsqrd_n=zeros(length(xm),1);
% rsqrd_a=zeros(length(xm),1);
% for i=1:length(xm);
%    mdl=fitlm(X(:,xm{i}),X(:,1));rsqrd_n(i)=mdl.Rsquared.Adjusted;
%    mdl=fitlm(X(:,xm{i}),X(:,2));rsqrd_a(i)=mdl.Rsquared.Adjusted;
% end
% 
% 
% % horizontal bar plot
% c0={'Forc','Vol','Sol','GHG'};
% c=categorical(c0,c0);
% figure
% barh(c,rsqrd_n)
% title('Forcing explain Nino3.4? (avg)')
% xlabel(['Adjusted r^2'],'fontsize',18)
% set(gca,'fontsize',16)
% 
% %print(['figs/linearfit_nino_forcings_r2.png'],'-dpng','-r300');
% 
% figure
% barh(c,rsqrd_a)
% title('Forcing explain AMO? (avg)')
% xlabel(['Adjusted r^2'],'fontsize',18)
% set(gca,'fontsize',16)
% 
% %print(['figs/linearfit_amo_forcings_r2.png'],'-dpng','-r300');




% %------------------------------------------
% % LOGISTIC REGRESSION ANALYSIS
% %------------------------------------------
% 
% % Bin the time series and do a logistric regression 
% 
% %ds=svmd; % only use megadroughts
% %ds=ds4; % use all decadal droughts
% ird_ann=[];
% for i=1:size(ds,1)
%   [~,~,irab] = intersect((ds(i,1)+analyrs1(1)-1):(ds(i,2)+analyrs1(1)-1),r_o:r_f); 
%   % Collect all years of drought
%   ird_ann=cat(1,ird_ann,irab);
% end
% 
% % Break up time series into single megadrought events and non-events
% 
% 
% y0=zeros(reconYrs,1);
% y0(ird_ann)=1;
% yd=y0(ir4);
% 
% % Break up time series by 'bl' year segments; find how many are drought years and bin accordingly
% %bl=10; % bin length; 21 = smallest bin length with no autocorrelation in the period 800-1600
% segs=fix(length(yd)/bl);
% fn=length(yd)-mod(length(yd),bl);
% ydr=reshape(yd(1:fn),bl,segs);
% ydi=find(sum(ydr)>=round(bl/2));
% yD=zeros(segs,1); % bins with not enough drought = 0
% yD(ydi)=1; % bins with more drought = 1
% 
% figure
% autocorr(yD)
% title('Is the drought time series autocorrelated?')
% 
% % Mean climate and forcing values in each bin
% X=zeros(segs,6);
% X0=xa_nino_ann(ir4)-mean(xa_nino_ann(ir2));
% X(1:segs,1)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,1)=min(reshape(X0(1:fn),bl,segs));
% 
% X0=xa_amo_mn(ir4)-mean(xa_amo_mn(ir2));
% X(1:segs,2)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,2)=sum(reshape(X0(1:fn),bl,segs));
% 
% X0=gmt_est_f(ig4)-mean(gmt_est_f(ig2));
% X(1:segs,3)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,3)=min(reshape(X0(1:fn),bl,segs));
% 
% X0=gmt_est_v(ig4)-mean(gmt_est_v(ig2));
% X(1:segs,4)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,4)=sum(reshape(X0(1:fn),bl,segs));
% 
% X0=gmt_est_s(ig4)-mean(gmt_est_s(ig2));
% X(1:segs,5)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,5)=sum(reshape(X0(1:fn),bl,segs));
% 
% X0=gmt_est_g(ig4)-mean(gmt_est_g(ig2));
% X(1:segs,6)=sum(reshape(X0(1:fn),bl,segs))./bl;
% %X(1:segs,6)=sum(reshape(X0(1:fn),bl,segs));
% 
% 
% Y=categorical(yD);
% %disp('logistic reg: Only using the ENSO and AMO as predictors')
% %[B,dev,stats] = mnrfit(X(:,1:2),Y);
% [B,dev,stats] = mnrfit(X,Y);
% 
% %close all
% 
% figure
% plot(1:length(stats.p),stats.p,'o')
% hold on
% plot(1:length(stats.p),0.05*ones(7,1))
% title('Logistic regression p-values')
% xticklabels({'Int.','Nino','AMO','Forc','Vol','Sol','GHG'});
% 
% %print(['figs/logregss_p_bsues.png'],'-dpng','-r300');
% 
% figure
% plot(stats.resid,'.')
% title('Logistic regression residuals')



% % Explore the predicted probabilities using the range of ENSO and AMO
% 
% % Define a grid of all possible Nino and AMO values
% %ns=xa_nino_ann(ir4,:)-nino_off;as=xa_amo_ens(ir4,:)-amo_off;
% % Just pick the mean values over which the model is built
% %ns=mean(xa_nino_ann(ir4,:)-nino_off,2);as=mean(xa_amo_ens(ir4,:)-amo_off,2);
% ns=X(:,1);as=X(:,2);
% [mx,my] = meshgrid(min(ns(:)):0.01:max(ns(:)),min(as(:)):0.01:max(as(:)));
% Xm(:,1)=mx(:);
% Xm(:,2)=my(:);
% 
% % Get probability for random observations
% pihat = mnrval(B,Xm);
% [pihat,dlow,dhi] = mnrval(B,Xm,stats);% get uncertainties
% % NOTE THAT UPPER AND LOWER BOUNDS ARE THE SAME!!
% 
% 
% % Reshape to match the grid
% pM=reshape(pihat(:,2),size(mx));
% pL=reshape(dlow(:,2),size(mx));
% %pU=reshape(dhi(:,2),size(mx));
% 
% 
% % Alternative plot with both probability and uncertainty
% figure
% hold on
% [C1,h1]=contour(mx,my,pM,'fill','on');
% set(gca,'tickdir','both')
% colormap(flip(gray))
% c=colorbar;
% c.Label.String = 'Est. decadal drought prob.';
% c.Label.FontSize = 16;
% [Cc,h]=contour(mx,my,pL,'showtext','on','linecolor',clrs.red,'linestyle','-.');
% clabel(Cc,h,'Color',clrs.red)
% set(gca,'fontsize',16)
% xlabel('Nino3.4 (10 yr anom. deg. C)','fontsize',18)
% ylabel('AMO (10 yr anom. deg. C)','fontsize',18)
% box on

%if prntfg=='y';print(['./figs/logreg_mdrought_prob_nino_amo_data&unc.png'],'-dpng','-r300');end
%print(['./figs/logreg_decdrought_prob_nino_amo_data&unc.png'],'-dpng','-r300');










% % FORCING AND DROUGHT 
% figure
% hold on
% %plot(analyrs3,DA_di(id3),'linewidth',1,'color',clrs.grey)
% %plot(analyrs2,Xa_di(ir2),'linewidth',1,'color',clrs.red)
% h1=plot(analyrs3,gmt_est_f(ig3),'linewidth',2,'color',clrs.black);
% h2=plot(analyrs3,gmt_est_v(ig3),'linewidth',2,'color',clrs.red);
% h3=plot(analyrs3,gmt_est_s(ig3),'linewidth',2,'color',clrs.blue);
% h4=plot(analyrs3,gmt_est_g(ig3),'linewidth',2,'color',clrs.green);
% %plot(analyrs2,smooth(Xa_di(ir2),10,'lowess'),'linewidth',2,'color',clrs.red)
% plot(analyrs3,zeros(length(analyrs3),1),'--','linewidth',2,'color',clrs.grey)
% set(gca,'fontsize',16)
% %ylim([-0.8 1])
% xlim([analyrs3(1) analyrs3(end)])
% xlabel('Year','fontsize',18)
% ylabel(['Forced response (anom. ' char(176) 'C)'],'fontsize',18)
% box on
% 
% 
% ax=gca;
% yMin=ax.YLim(1)+0.5;yMax=ax.YLim(2)-0.5;
% %y=gmt_est_mn(ig3); % using 'l' from above
% %yMin=min(y)-0.5;yMax=max(y)+0.5;
% ylim([yMin yMax])
% % Add a patches for megadroughts
% for i=1:length(ds)
%     p=patch(analyrs3(1)+[ds(i,1) ds(i,2) ds(i,2) ds(i,1)],[yMin yMin yMax yMax],clrs.orange);
%     p.FaceAlpha=0.35;p.EdgeColor='none';
% end
% 
% % Add patches for all other decadal droughts
% % for i=10:size(ds,1)
% %     p=patch(analyrs3(1)+[ds(i,1) ds(i,2) ds(i,2) ds(i,1)],[yMin yMin yMax yMax],clrs.olive);
% %     p.FaceAlpha=0.25;p.EdgeColor='none';
% % end
% 
% legend([h1,h2,h3,h4],{'All','Volc.','Solar','GHG'},'location','southwest')
% hold off
% 
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]
% 
% %print(['figs/forcings_alldroughts_800_2000.png'],'-dpng','-r300');
% 
% % FORCING, NINO, AMO, AND DROUGHT 
% figure
% hold on
% %plot(analyrs3,DA_di(id3),'linewidth',1,'color',clrs.grey)
% %plot(analyrs2,Xa_di(ir2),'linewidth',1,'color',clrs.red)
% h1=plot(mon_r(irm3_a:irm3_b),xa_nino_mn(irm3_a:irm3_b)-nino_off,'linewidth',1.5,'color',clrs.blue);
% h3=plot(analyrs3,xa_amo_mn(ir3)-amo_off,'linewidth',2,'color',clrs.green);
% h2=plot(analyrs3,gmt_est_f(ig3),'linewidth',2,'color',clrs.red);
% plot(analyrs3,zeros(length(analyrs3),1),'--','linewidth',2,'color',clrs.yellow)
% set(gca,'fontsize',12)
% %ylim([-0.8 1])
% xlim([analyrs3(1) analyrs3(end)])
% xlabel('Year','fontsize',13)
% ylabel(['Temperature anomaly (' char(176) 'C)'],'fontsize',13)
% box on
% 
% ax=gca;
% yMin=ax.YLim(1)+0.5;yMax=ax.YLim(2)-0.5;
% %y=gmt_est_mn(ig3); % using 'l' from above
% %yMin=min(y)-0.5;yMax=max(y)+0.5;
% ylim([yMin yMax])
% % Add a patches for megadroughts
% for i=1:length(ds)
%     p=patch(analyrs3(1)+[ds(i,1) ds(i,2) ds(i,2) ds(i,1)],[yMin yMin yMax yMax],clrs.orange);
%     p.FaceAlpha=0.35;p.EdgeColor='none';
% end
% %legend([h1,h3,h2],{'Nino3.4','AMO','Forcing'},'location','southeast')
% hold off
% 
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 8 3]; % [left bottom width height]
% 
% %print(['figs/mdrought_nino_forcing_amo_800_2000.png'],'-dpng','-r300');

















% % Bivariate logistic regression with interaction
% %Xi(:,1)=mean(xa_nino_ann(ir4,:)-nino_off,2);
% %Xi(:,2)=mean(xa_amo_ens(ir4,:)-amo_off,2);
% %Xi(:,3)=Xi(:,1).*Xi(:,2);
% %[Bi,devi,statsi] = mnrfit(Xi,Y);
% 
% 
% 
% % Explore the predicted probabilities using the range of ENSO and AMO
% 
% % Define a grid of all possible Nino and AMO values
% %ns=xa_nino_ann(ir4,:)-nino_off;as=xa_amo_ens(ir4,:)-amo_off;
% % Just pick the mean values over which the model is built
% ns=mean(xa_nino_ann(ir4,:)-nino_off,2);as=mean(xa_amo_ens(ir4,:)-amo_off,2);
% [mx,my] = meshgrid(min(ns(:)):0.01:max(ns(:)),min(as(:)):0.01:max(as(:)));
% Xm(:,1)=mx(:);
% Xm(:,2)=my(:);
% 
% % Get probability for random observations
% pihat = mnrval(B,Xm);
% [pihat,dlow,dhi] = mnrval(B,Xm,stats);% get uncertainties
% % NOTE THAT UPPER AND LOWER BOUNDS ARE THE SAME!!
% 
% 
% % Reshape to match the grid
% pM=reshape(pihat(:,2),size(mx));
% pL=reshape(dlow(:,2),size(mx));
% %pU=reshape(dhi(:,2),size(mx));
% 
% % figure
% % [Cc,h]=contour(mx,my,pM,'ShowText','on','fill','on','linecolor','r');
% % colormap(flip(gray))
% % clabel(Cc,h,'FontSize',14,'Color','red')
% % title('Predicted megadrought probability','fontsize',20)
% % set(gca,'fontsize',16)
% % xlabel('Nino3.4 (annual anom. deg. C)','fontsize',18)
% % ylabel('AMO (anom. deg. C)','fontsize',18)
% % 
% % if prntfg=='y';print(['./figs/logreg_mdrought_prob_nino_amo.png'],'-dpng','-r300');end
% % 
% % 
% % figure
% % [Cc,h]=contour(mx,my,pL,'ShowText','on','linecolor','k');
% % clabel(Cc,h,'FontSize',14,'Color','black')
% % title('Probability 95% conf.','fontsize',20)
% % set(gca,'fontsize',16)
% % xlabel('Nino3.4 (annual anom. deg. C)','fontsize',18)
% % ylabel('AMO (anom. deg. C)','fontsize',18)
% % 
% % if prntfg=='y';print(['./figs/logreg_mdrought_prob_nino_amo_unc.png'],'-dpng','-r300');end
% 
% 
% % Alternative plot with both probability and uncertainty
% figure
% hold on
% [C1,h1]=contour(mx,my,pM,'fill','on');
% set(gca,'tickdir','both')
% colormap(flip(gray))
% c=colorbar;
% c.Label.String = 'Est. megadrought probability';
% c.Label.FontSize = 16;
% [Cc,h]=contour(mx,my,pL,'showtext','on','linecolor',clrs.red,'linestyle','-.');
% clabel(Cc,h,'Color',clrs.red)
% set(gca,'fontsize',16)
% xlabel('Nino3.4 (annual anom. deg. C)','fontsize',18)
% ylabel('AMO (anom. deg. C)','fontsize',18)
% box on
% 
% %if prntfg=='y';print(['./figs/logreg_mdrought_prob_nino_amo_data&unc.png'],'-dpng','-r300');end
% 
















% % AMO offset
% amo_off=mean(mean(xa_amo_ens(ir2,:),2));
% xa_amo_mn=mean(xa_amo_ens,2);
% 
% % FORCING AND DROUGHT 
% figure
% hold on
% %plot(analyrs3,DA_di(id3),'linewidth',1,'color',clrs.grey)
% %plot(analyrs2,Xa_di(ir2),'linewidth',1,'color',clrs.red)
% h1=plot(analyrs3,xa_amo_mn(ir3)-amo_off,'linewidth',2,'color',[0 0 0]);
% plot(analyrs3,gmt_est_mn(ig3),'linewidth',2,'color',clrs.red)
% %plot(analyrs2,smooth(Xa_di(ir2),10,'lowess'),'linewidth',2,'color',clrs.red)
% plot(analyrs3,zeros(length(analyrs3),1),'--','linewidth',2,'color',clrs.yellow)
% set(gca,'fontsize',16)
% %ylim([-0.8 1])
% xlim([analyrs3(1) analyrs3(end)])
% xlabel('Year','fontsize',18)
% ylabel(['Temperature anomaly (' char(176) 'C)'],'fontsize',18)
% box on
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]
% 
% y=Xa_di(ir3); % using 'l' from above
% yMin=min(y)-0.5;yMax=max(y)+0.5;
% ylim([yMin yMax])
% % Add a patches for top five times from truth
% for i=1:9
%     p=patch(analyrs3(1)+[ds4(i,1) ds4(i,2) ds4(i,2) ds4(i,1)],[yMin yMin yMax yMax],clrs.blue);
%     p.FaceAlpha=0.35;p.EdgeColor='none';
% end
% 
% %legend({'NADA','NADA-11'},'location','southeast')
% hold off
% 
% %print(['figs/gmtest_drght_AMO_cmp.png'],'-dpng','-r300');
% 
% 
% 
% 
% 
% nino_off=mean(mean(xa_nino_ens(irm2_a:irm2_b,:),2));
% 
% 
% % FORCING AND DROUGHT 
% figure
% hold on
% %plot(analyrs3,DA_di(id3),'linewidth',1,'color',clrs.grey)
% %plot(analyrs2,Xa_di(ir2),'linewidth',1,'color',clrs.red)
% h1=plot(mon_r(irm3_a:irm3_b),xa_nino_mn(irm3_a:irm3_b)-nino_off,'linewidth',2,'color',[0 0 0]);
% plot(analyrs3,gmt_est_mn(ig3),'linewidth',2,'color',clrs.red)
% %plot(analyrs2,smooth(Xa_di(ir2),10,'lowess'),'linewidth',2,'color',clrs.red)
% plot(analyrs3,zeros(length(analyrs3),1),'--','linewidth',2,'color',clrs.yellow)
% set(gca,'fontsize',16)
% %ylim([-0.8 1])
% xlim([analyrs3(1) analyrs3(end)])
% xlabel('Year','fontsize',18)
% ylabel(['Temperature anomaly (' char(176) 'C)'],'fontsize',18)
% box on
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]
% 
% y=Xa_di(ir3); % using 'l' from above
% yMin=min(y)-0.5;yMax=max(y)+0.5;
% ylim([yMin yMax])
% % Add a patches for top five times from truth
% for i=1:9
%     p=patch(analyrs3(1)+[ds4(i,1) ds4(i,2) ds4(i,2) ds4(i,1)],[yMin yMin yMax yMax],clrs.blue);
%     p.FaceAlpha=0.35;p.EdgeColor='none';
% end
% 
% %legend({'NADA','NADA-11'},'location','southeast')
% hold off
% 
% %print(['figs/gmtest_drght_Nino_cmp.png'],'-dpng','-r300');
% 
% 









