%Caculates principal components and averages from the NOAMER
%data-set using the normalizations discussed in the comment on MM05.

clear;

%Load North American Tree Rings
NOAMER;
proxy_time=proxy(:,2);
proxy=proxy(:,3:end);
pnum=length(proxy(1,:));

%Load Northern Hemisphere instrumental temperatures
Northern_Hemisphere_Temperature;
NH_time=NH_temp(1:2:end,1);
NH_temp=NH_temp(1:2:end,end);

%Define empty matrices
pMBH=[];
pfull=[];
pMM=[];

%Normalize the tree rings
for ct=1:pnum,
  x=proxy(:,ct);

  %The MM05 normalization
  pMM(:,end+1)=x-mean(x);
    
  %The MBH98 normalization
  pl=find(proxy_time>=1902 & proxy_time<=1980 & isnan(x)==0);
  x=(x-mean(x(pl)))/std(detrend(x(pl))); 
  %x=(x-mean(x(pl)))/std(x(pl)); %Uncoment to see result without detrending.
  pMBH(:,end+1)=x; 
  
  %Full normalization
  pl=find(proxy_time>=1400 & proxy_time<=1980 & isnan(x)==0);
  x=(x-mean(x(pl)))/std(x(pl)); 
  pfull(:,end+1)=x; 
end;

%Simple averages
avg(1,:)=mean(pMBH');
avg(2,:)=mean(pMM');
avg(3,:)=mean(pfull');

%PCA and fraction of variance explained
[U S pc1(1,:)]=svds(pMBH',1);   
[U S pc1(2,:)]=svds(pMM',1);    
[U S pc1(3,:)]=svds(pfull',1);  

[U S pc2]=svds(pMBH',3);    Spc(1,:)=sum(S).^2/sum(sum(pMBH.^2));
[U S pc2]=svds(pMM',3);     Spc(2,:)=sum(S).^2/sum(sum(pMM.^2));
[U S pc2]=svds(pfull',3);   Spc(3,:)=sum(S).^2/sum(sum(pfull.^2));


%Orient pc1 according to averages
for ct=1:3,
  pc1(ct,:)=pc1(ct,:)*sign(xcPH(pc1(ct,:),avg(ct,:),1));
end;


%-----------------------------------Variances (Figure 1)

figure(1); clf; subplot(5,1,[1 2 3]); hold on;
col=[0.75 0.75 0.75];

psMBH=sum(pMBH.^2)/sum(sum(pMBH.^2));
[dummy j]=sort(psMBH);
plot(psMBH(j),'k+');
psfull=sum(pfull.^2)/sum(sum(pfull.^2));
plot(psfull(j),'k.');
psMM=sum(pMM.^2)/sum(sum(pMM.^2));
plot(psMM(j),'kx');
h=axis; 

%The below indices correspond with records listed in Table 1 of MM05.
recs=[7 13 14 15 18 19 27 28 29 30 51 53 54 55 57];
for ct=1:length(recs),
 pl=find(j==recs(ct));
 h2=fill([pl-.5 pl+.5 pl+.5 pl-.5],[h(3) h(3) h(4) h(4)],col);
 set(h2,'edgecolor',col);
end;

plot(psMBH(j),'k+');
psfull=sum(pfull.^2)/sum(sum(pfull.^2));
plot(psfull(j),'k.');
psMM=sum(pMM.^2)/sum(sum(pMM.^2));
plot(psMM(j),'kx');
h=xlabel('record number'); 
h=ylabel('fraction of total variance'); 

axis tight;
print -depsc ../Reply/Figures/Variance.eps;


 
%-------------------Compute various statistics

%variance ratios

recs_other=1:70;
for ct=1:length(recs),
  pl=find(recs_other~=recs(ct));
  recs_other=recs_other(pl);
end;
ratioMBH=mean(psMBH(recs))./mean(psMBH(recs_other)),
ratioMM=mean(psMM(recs))./mean(psMM(recs_other)),

%hockey stick index
pl=find(proxy_time>=1902); 
for ct=1:3,
  HSI(ct)=(mean(pc1(ct,pl))-mean(pc1(ct,:)))/std(pc1(ct,:));
end;

%cross-correlations between averages and PC1s 
for ct1=1:3,
  for ct2=1:3,
    xc_avg_pc1(ct1,ct2)=xcPH(avg(ct1,:),pc1(ct2,:),1);
  end;
end;



%--------------------------------Averages and PCA (Figure 2)


%Smooth records using a 11 year Hamming window
sm_year=11;
for ct=1:3,
  avg(ct,:)=smoothPH(avg(ct,:),sm_year);
  pc1(ct,:)=smoothPH(pc1(ct,:),sm_year);
end;
NH_temp=smoothPH(NH_temp,sm_year);

%Normalize averages and principal components to the instrumental record
index_instr=find(NH_time>=1902 & NH_time<=1980);
mean_inst=mean(NH_temp(index_instr));
std_inst=std(NH_temp(index_instr));
index_proxy=find(proxy_time>=1902 & proxy_time<=1980);
for ct=1:3,
  avg(ct,:)=(avg(ct,:)-mean(avg(ct,index_proxy)))*std_inst/std(avg(ct,index_proxy))+mean_inst;
  pc1(ct,:)=(pc1(ct,:)-mean(pc1(ct,index_proxy)))*std_inst/std(pc1(ct,index_proxy))+mean_inst;
end;

%Make figure of averages and pc1s
labels='abc';
figure(2); clf;
for ct=1:3,
  subplot(3,1,ct); hold on;
  plot([1400,2000],[0 0],'k--');
  h=plot(NH_time,NH_temp,'k-.'); 
  set(h,'linewidth',1.5);
  h=plot(proxy_time,avg(ct,:),'k'); 
  set(h,'color',col,'linewidth',1.5);
  plot(proxy_time,pc1(ct,:),'k');
  axis tight; 
  h=text(2000,0.25,['(',labels(ct),')']);
  h=ylabel('anomalies (^\circC)'); 
end;

%orient tall;
%print -depsc ../Reply/Figures/Reconstructions.eps;