
% 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=800:1600;
analyrs2=1601:1925;
analyrs3=800:2000;
analyrs4=800:1925;
[~,~,ir1] = intersect(analyrs1,r_o:r_f); % indices of recon segment
[~,~,ir2] = intersect(analyrs2,r_o:r_f); % indices of recon segment
[~,~,ir3] = intersect(analyrs3,r_o:r_f); % indices of recon segment
[~,~,ir4] = intersect(analyrs4,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;

ovrlprng=rng_intr([r_o,r_f],[analyrs2(1),analyrs2(end)]);
mon_r=r_o:(1/12):(r_f+11/12);
irm2_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
irm2_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;

ovrlprng=rng_intr([r_o,r_f],[analyrs3(1),analyrs3(end)]);
mon_r=r_o:(1/12):(r_f+11/12);
irm3_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
irm3_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(ir2,:),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(irm2_a:irm2_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);


% PULL OUT SPATIAL MEAN OVER AMERICAN SOUTHWEST
%load('/d1/nsteiger/climate-data/cesm-lme/Nrth_A_Swthwst_mask.mat')
varpth=['/d6/haibo/lme/cam5/BLMTRC5CN/ALLFORCEINGS/010'];
sftlf=permute(ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'landmask'),[2 1]);
%lat=ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'lat');
%lon=ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'lon');
% NASW = box of 235:255 E (125:105 W); 31:42 N (just the US, not Mexico)
% Find nearest lat/lon equal to or outside bounds
lonext=[235 255];latext=[31 42];
ln_a=find((abs(lon-lonext(1))-min(abs(lon-lonext(1))))<1e-5);
ln_b=find((abs(lon-lonext(2))-min(abs(lon-lonext(2))))<1e-5);
lt_a=find((abs(lat-latext(1))-min(abs(lat-latext(1))))<1e-5);
lt_b=find((abs(lat-latext(2))-min(abs(lat-latext(2))))<1e-5);
nasw=zeros(size(sftlf));
nasw(lt_a:lt_b,ln_a:ln_b)=sftlf(lt_a:lt_b,ln_a:ln_b);
nasw(nasw==0)=NaN;

% 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,nasw);
t2m_nasw=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')




%======================
%======  PHYDA  =======
%======================

% k=find(state_tp=='q');
% lat=xmeta{k}.lat;
% lon=xmeta{k}.lon;
% 
% % PULL OUT SPATIAL MEAN OVER AMERICAN SOUTHWEST
% %load('/d1/nsteiger/climate-data/cesm-lme/Nrth_A_Swthwst_mask.mat')
% varpth=['/d6/haibo/lme/cam5/BLMTRC5CN/ALLFORCEINGS/010'];
% sftlf=permute(ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'landmask'),[2 1]);
% %lat=ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'lat');
% %lon=ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'lon');
% % NASW = box of 235:255 E (125:105 W); 31:42 N (just the US, not Mexico)
% % Find nearest lat/lon equal to or outside bounds
% lonext=[235 255];latext=[31 42];
% ln_a=find((abs(lon-lonext(1))-min(abs(lon-lonext(1))))<1e-5);
% ln_b=find((abs(lon-lonext(2))-min(abs(lon-lonext(2))))<1e-5);
% lt_a=find((abs(lat-latext(1))-min(abs(lat-latext(1))))<1e-5);
% lt_b=find((abs(lat-latext(2))-min(abs(lat-latext(2))))<1e-5);
% nasw=zeros(size(sftlf));
% nasw(lt_a:lt_b,ln_a:ln_b)=sftlf(lt_a:lt_b,ln_a:ln_b);
% nasw(nasw==0)=NaN;
% 
% % Latitude weighting
% A=cosd(repmat(xmeta{k}.lat,[1 length(xmeta{k}.lon)]));
% 
% % PDSI ensemble
% Xa_di_ens=zeros(reconYrs,sub_ens);
% for i=1:sub_ens
%    Xvar=squeeze(Xa_ens(xmeta{k}.id_X(1):xmeta{k}.id_X(2),i,:));
%    Xpdsi=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),reconYrs);
%    Xpdsi(xmeta{k}.lndidx,:)=Xvar;
%    Xpdsi=reshape(Xpdsi,length(xmeta{k}.lat),length(xmeta{k}.lon),reconYrs);
%    Xv_msk=bsxfun(@times,Xpdsi,nasw);
%    Xa_di_ens(:,i)=wmean_a(Xv_msk,A); % RECONSTRUCTION DROUGHT INDEX
% end
% 
% % Get the ensemble mean; identical to the original reconstruction
% Xvar=squeeze(Xa_m(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:));
% Xpdsi=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),reconYrs);
% Xpdsi(xmeta{k}.lndidx,:)=Xvar;
% Xpdsi=reshape(Xpdsi,length(xmeta{k}.lat),length(xmeta{k}.lon),reconYrs);
% Xv_msk=bsxfun(@times,Xpdsi,nasw);
% Xa_di=wmean_a(Xv_msk,A); % RECONSTRUCTION DROUGHT INDEX
% 
% % How to compute the megadroughts? 1 => 11 year running mean, 2 => 2S2E
% %ds1=mdrought(Xa_di(ir1),1,'y'); % 'y' indicates that the mean is removed
% %ds2=mdrought(Xa_di(ir1),2,'y');
% 
% % How to compute the megadroughts? 1 => 11 year running mean, 2 => 2S2E
% ds4=mdrought(Xa_di(ir4),1,'y'); % 'y' indicates that the mean is removed
% % Most severe droughts have severity < -10 and all occur before 1600 (last one ends in 1590)
% si=find(ds4(:,3)<-10);
% %si=find(ds4(:,2)+analyrs1(1)<analyrs1(end));
% svmd=ds4(si,:);


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

load('svmd_ens_90pagree.mat')
%load('svmd_ens_50pagree.mat')
ds=svmd;

%======================
%======  NADA  ========
%======================

% da=permute(ncread('/d1/nsteiger/drought-atlas/nada_hd2_cl.nc','pdsi'),[2 3 1]);
% %nada(nada<-99)=NaN; % put nans in fill spots
% lat_da=ncread('/d1/nsteiger/drought-atlas/nada_hd2_cl.nc','lat');
% lon_da=ncread('/d1/nsteiger/drought-atlas/nada_hd2_cl.nc','lon');
% yrs_da=ncread('/d1/nsteiger/drought-atlas/nada_hd2_cl.nc','time');
% 
% 
% % Coarsen NADA to match climate model grid
% ndm=load('nada_mask.mat');
% % create grids for interpolation
% nlns=find(lon_da<0);lon_da(nlns)=lon_da(nlns)+360;
% [xn,yn]=meshgrid(lon_da,lat_da);
% [xc,yc]=meshgrid(lon,lat);
% % note masks of both data types
% nm=find(~isnan(da(:,:,end)));
% cm=find(ndm.sftlf==1);
% da(da<-99)=NaN; % put nans in fill spots but after determining mask
% 
% da_i=nan(length(lat),length(lon),size(da,3));
% for i=1:size(da,3)
%    % linear interpolation and extrapolation for irregular data
%    da1=da(:,:,i);
%    F = scatteredInterpolant(xn(nm),yn(nm),da1(nm),'linear','linear');
%    vq=F(xc(cm),yc(cm)); % do the interpolation
%    % expand out the interpolated data points onto the climate model mask
%    vd=nan(length(lat),length(lon));
%    vd(cm)=vq; % full interpolated field
%    da_i(:,:,i)=vd;
% end
% 
% 
% 
% % PULL OUT SPATIAL MEAN OVER AMERICAN SOUTHWEST
% %load('/d1/nsteiger/climate-data/cesm-lme/Nrth_A_Swthwst_mask.mat')
% varpth=['/d6/haibo/lme/cam5/BLMTRC5CN/ALLFORCEINGS/010'];
% sftlf=permute(ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'landmask'),[2 1]);
% %lat=ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'lat');
% %lon=ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'lon');
% % NASW = box of 235:255 E (125:105 W); 31:42 N (just the US, not Mexico)
% % Find nearest lat/lon equal to or outside bounds
% lonext=[235 255];latext=[31 42];
% ln_a=find((abs(lon-lonext(1))-min(abs(lon-lonext(1))))<1e-5);
% ln_b=find((abs(lon-lonext(2))-min(abs(lon-lonext(2))))<1e-5);
% lt_a=find((abs(lat-latext(1))-min(abs(lat-latext(1))))<1e-5);
% lt_b=find((abs(lat-latext(2))-min(abs(lat-latext(2))))<1e-5);
% nasw=zeros(size(sftlf));
% nasw(lt_a:lt_b,ln_a:ln_b)=sftlf(lt_a:lt_b,ln_a:ln_b);
% 
% nasw(nasw==0)=NaN;
% A=cosd(repmat(xmeta{k}.lat,[1 length(xmeta{k}.lon)]));
% 
% DA_msk=bsxfun(@times,da_i,nasw);
% DA_di=wmean_a(DA_msk,A); % DROUGHT ATLAS DROUGHT INDEX
% % Indices of the drought atlas
% analyrs5=1870:2000;
% [~,~,id1] = intersect(analyrs1,yrs_da); % indices of recon segment
% [~,~,id2] = intersect(analyrs2,yrs_da); % indices of recon segment
% [~,~,id3] = intersect(analyrs3,yrs_da); % indices of recon segment
% [~,~,id4] = intersect(analyrs4,yrs_da); % indices of recon segment
% [~,~,id5] = intersect(analyrs5,yrs_da); % indices of recon segment
% 
% [~,~,ir5] = intersect(analyrs5,r_o:r_f); % indices of recon segment
% 
% 
% % MEGADROUGHTS IN NADA	
% nada_ds4=mdrought(DA_di(ir4),1,'y'); % 'y' indicates that the mean is removed



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

% Full forcing
load('../megad-forcing/output_twoboxCM_full_06-Aug-2018_15:09:13.mat','gmt_fr','myrs')

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

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

% solar 
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')
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')
gmt_est_mn=mean(gmt_fr,2);
gmt_est_g=gmt_est_mn-mean(gmt_est_mn(ig2));

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

load distinct_colors_20.mat

% Analysis options
plot_p='y'; % plot proxies?
prntfg='n'; % Save figures?



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

% 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'};
c=categorical(c0,c0);


% Averaging length
%al=10;

bs=50; disp('Reduced bootstrap iterations...')
%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),6,bs,length(al));
p_bs=NaN(length(xm),6,bs,length(al));

k=1;
for j=al 
   % estimate NASW PDSI
   clear X4 y4
   y4=detrend(ann2xt(Xa_di(ir4),j));
   X4(:,1)=detrend(ann2xt(xa_nino_ann(ir4),j));
   X4(:,2)=detrend(ann2xt(xa_amo_mn(ir4),j));
   X4(:,3)=detrend(ann2xt(gmt_est_f(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=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),6,bs*length(al));
p0=reshape(p_bs,length(xm),6,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,:));



% Distributions of the r2 values
figure
hold on
boxplot(r2',c)
plot(1:length(c),zeros(length(c),1),'--','color',clrs.grey,'linewidth',2)
ylabel('Adjusted r^2','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/linear_r2_jackknife.png'],'-dpng','-r300');
%print(['figs/linear_r2_bootstrap.png'],'-dpng','-r300');


% Model coefficients
figure
hold on
boxplot(b1,c(1:5))
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 4 4]; % [left bottom width height]

%print(['figs/linear_betas_jackknife.png'],'-dpng','-r300');
%print(['figs/linear_betas_bootstrap.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','northeast')
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.png'],'-dpng','-r300');

close all

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


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

% Averaging length
%al=10;

bs=500;
%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(ir4),j));
   y4a0=detrend(ann2xt(xa_amo_mn(ir4),j));
   y4t0=detrend(ann2xt(t2m_nasw(ir4),j));
   X4(:,1)=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);
   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,:));


% 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(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)
ylim([-0.65 0.8])
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_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)
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');




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









