%%%% Matlab version of Hockey Stick 
%
%Mann, M. E., R. S. Bradley, and M. K. Hughes (1998), Global-scale
%temperature patterns and climate forcing over the past six centuries,
%Nature, 392, 779– 787.
%
%Mann, M. E., R. S. Bradley, and M. K. Hughes (1999), Northern
%Hemisphere temperatures during the past millennium: Inferences, uncertainties,
%and limitations, Geophys. Res. Lett., 26, 759– 762.
%
% Tested with Matlab version 7.4. (R2007a)
% UC - uc_edit at yahoo.com  - http://signals.auditblogs.com/
%%%%% hockeystick.txt (.m) ver 1.1 July 08


close all
%%% PART I       %%%  
%%% Get the data from web and write to files (once)%%%
%%% Your firewall might have a word at this point..%

if ~exist('downloaded_all.mat','file')

% Original reconstruction
urlwrite('http://holocene.meteo.psu.edu/shared/research/ONLINE-PREPRINTS/Millennium/DATA/RECONS/nhem-recon.dat','nhem_recon.dat')

% MBH98 rec and Instrumental mean series

urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/FigureData/nhmean.txt','nhmean.txt')

%  
%urlwrite('ftp://ftp.ncdc.noaa.gov/pub/data/paleo/paleocean/by_contributor/mann1998/nhem-dense.dat','nhem_dense.dat');

% Sparse instrumental data

urlwrite('http://www.ncdc.noaa.gov/paleo/ei/ei_data/nhem-sparse.dat','nhem_sparse.dat')

% gridpoints

urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/gridpoints.txt','gridpoints.txt')

% proxies

urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1400.txt','data1400.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1450.txt','data1450.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1500.txt','data1500.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1600.txt','data1600.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1700.txt','data1700.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1730.txt','data1730.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1750.txt','data1750.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1760.txt','data1760.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1780.txt','data1780.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1800.txt','data1800.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1820.txt','data1820.txt')


% MBH99 AD1000

urlwrite('http://holocene.meteo.psu.edu/shared/research/ONLINE-PREPRINTS/Millennium/DATA/PROXIES/pc1-fixed.dat','pc1_fixed.dat');

% http://holocene.meteo.psu.edu/shared/research/ONLINE-PREPRINTS/Millennium/DATA/PROXIES/
% morc014 is missing, so using CA version

urlwrite('http://www.climateaudit.org/data/mbh99/proxy.txt','proxy.txt')


% Eigenvalues for S-matrix
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/tpca-eigenvals.txt','tpca_eigenvals.txt')



% U matrix

urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc01.txt','pc01.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc02.txt','pc02.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc03.txt','pc03.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc04.txt','pc04.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc05.txt','pc05.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc06.txt','pc06.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc07.txt','pc07.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc08.txt','pc08.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc09.txt','pc09.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc10.txt','pc10.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc11.txt','pc11.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc12.txt','pc12.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc13.txt','pc13.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc14.txt','pc14.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc15.txt','pc15.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc16.txt','pc16.txt')

% V matrix

urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof01.txt','eof01.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof02.txt','eof02.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof03.txt','eof03.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof04.txt','eof04.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof05.txt','eof05.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof06.txt','eof06.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof07.txt','eof07.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof08.txt','eof08.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof09.txt','eof09.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof10.txt','eof10.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof11.txt','eof11.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof12.txt','eof12.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof13.txt','eof13.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof14.txt','eof14.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof15.txt','eof15.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof16.txt','eof16.txt')

% Standard deviations

urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/tpca-standard.txt','tpca_standard.txt')


downloaded_all=1;

save downloaded_all downloaded_all % just to mark successful download


end % if exist



%%% PART II       %%%  
%%% Load the data from files and preprocess for calibration%%%


load nhem_recon.dat

% Proxies

fid1=fopen('proxy.txt','r');
dds=[];
tline=fgetl(fid1);
while 1
            tline = fgetl(fid1);
            if ~ischar(tline), break, end
            
            if tline(1:3)=='976', break,end % NAs
            
            dds=[dds;str2num(tline)];
            
end

temp=str2num(tline(1:14)); temp2=str2num(tline(19:end)); % replace NAs with previous values (kind of manually..) 
dds=[dds; temp 1.063 temp2];

tline=fgetl(fid1);
temp=str2num(tline(1:13)); temp2=str2num(tline(18:end)); % replace NAs with previous values
dds=[dds; temp 1.063 temp2];

tline=fgetl(fid1);
temp=str2num(tline(1:14)); temp2=str2num(tline(19:end)); % replace NAs with previous values
dds=[dds; temp 1.063 temp2];

tline=fgetl(fid1);
temp=str2num(tline(1:14)); temp2=str2num(tline(19:end)); % replace NAs with previous values
dds=[dds; temp 1.063 temp2];

tline=fgetl(fid1);
temp=str2num(tline(1:13)); temp2=str2num(tline(18:end)); % replace NAs with previous values
dds=[dds; temp 1.063 temp2];

tline=fgetl(fid1);
temp=str2num(tline(1:13)); temp2=str2num(tline(18:end)); % replace NAs with previous values
dds=[dds; temp 1.063 temp2];



fclose(fid1);

ad1000=(dds(:,3:16));


ad1000b=ad1000; % fixed PC1 version
load pc1_fixed.dat
ad1000b(:,3)=pc1_fixed(:,2);

% standardize by detrended std  ( http://www.nature.com/nature/journal/v430/n6995/extref/METHODS/README.txt )   

AD1000=(ad1000-repmat(mean(ad1000(903:end,:)),981,1))./repmat(std(detrend(ad1000(903:end,:))),981,1);
AD1000b=(ad1000b-repmat(mean(ad1000b(903:end,:)),981,1))./repmat(std(detrend(ad1000b(903:end,:))),981,1);


load data1400.txt


ad1400=[zeros(400,22); data1400(:,2:end)];  % fix size to 981, helps later

for ii=1:981
    iNaN=isnan(ad1400(ii,:));  %% if there are NaNs, replace with previous 
    if sum(iNaN)
        fiNaN=find(iNaN);
        ad1400(ii,fiNaN)=ad1400(ii-1,fiNaN);
    end
end
    
AD1400=(ad1400-repmat(mean(ad1400(903:end,:)),981,1))./repmat(std(detrend(ad1400(903:end,:))),981,1);


load data1450.txt
ad1450=[zeros(450,25); data1450(:,2:end)];  

for ii=1:981
    iNaN=isnan(ad1450(ii,:));  %% if there are NaNs, replace with previous 
    if sum(iNaN)
        fiNaN=find(iNaN);
        ad1450(ii,fiNaN)=ad1450(ii-1,fiNaN);
    end
end


AD1450=(ad1450-repmat(mean(ad1450(903:end,:)),981,1))./repmat(std(detrend(ad1450(903:end,:))),981,1);


load data1500.txt
ad1500=[zeros(500,28); data1500(:,2:end)];  

for ii=1:981
    iNaN=isnan(ad1500(ii,:));  %% if there are NaNs, replace with previous 
    if sum(iNaN)
        fiNaN=find(iNaN);
        ad1500(ii,fiNaN)=ad1500(ii-1,fiNaN);
    end
end


AD1500=(ad1500-repmat(mean(ad1500(903:end,:)),981,1))./repmat(std(detrend(ad1500(903:end,:))),981,1);

load data1600.txt
ad1600=[zeros(600,57); data1600(:,2:end)];

for ii=1:981
    iNaN=isnan(ad1600(ii,:));  %% if there are NaNs, replace with previous 
    if sum(iNaN)
        fiNaN=find(iNaN);
        ad1600(ii,fiNaN)=ad1600(ii-1,fiNaN);
    end
end


AD1600=(ad1600-repmat(mean(ad1600(903:end,:)),981,1))./repmat(std(detrend(ad1600(903:end,:))),981,1);

load data1700.txt
ad1700=[zeros(700,74); data1700(:,2:end)];

for ii=1:981
    iNaN=isnan(ad1700(ii,:));  %% if there are NaNs, replace with previous 
    if sum(iNaN)
        fiNaN=find(iNaN);
        ad1700(ii,fiNaN)=ad1700(ii-1,fiNaN);
    end
end


AD1700=(ad1700-repmat(mean(ad1700(903:end,:)),981,1))./repmat(std(detrend(ad1700(903:end,:))),981,1);


load data1730.txt
ad1730=[zeros(730,79); data1730(:,2:end)];

for ii=1:981
    iNaN=isnan(ad1730(ii,:));  %% if there are NaNs, replace with previous 
    if sum(iNaN)
        fiNaN=find(iNaN);
        ad1730(ii,fiNaN)=ad1730(ii-1,fiNaN);
    end
end

AD1730=(ad1730-repmat(mean(ad1730(903:end,:)),981,1))./repmat(std(detrend(ad1730(903:end,:))),981,1);

load data1750.txt
ad1750=[zeros(750,89); data1750(:,2:end)];

for ii=1:981
    iNaN=isnan(ad1750(ii,:));  %% if there are NaNs, replace with previous 
    if sum(iNaN)
        fiNaN=find(iNaN);
        ad1750(ii,fiNaN)=ad1750(ii-1,fiNaN);
    end
end

AD1750=(ad1750-repmat(mean(ad1750(903:end,:)),981,1))./repmat(std(detrend(ad1750(903:end,:))),981,1);

load data1760.txt
ad1760=[zeros(760,93); data1760(:,2:end)];
for ii=1:981
    iNaN=isnan(ad1760(ii,:));  %% if there are NaNs, replace with previous 
    if sum(iNaN)
        fiNaN=find(iNaN);
        ad1760(ii,fiNaN)=ad1760(ii-1,fiNaN);
    end
end

AD1760=(ad1760-repmat(mean(ad1760(903:end,:)),981,1))./repmat(std(detrend(ad1760(903:end,:))),981,1);

load data1780.txt
ad1780=[zeros(780,97); data1780(:,2:end)];

for ii=1:981
    iNaN=isnan(ad1780(ii,:));  %% if there are NaNs, replace with previous 
    if sum(iNaN)
        fiNaN=find(iNaN);
        ad1780(ii,fiNaN)=ad1780(ii-1,fiNaN);
    end
end

AD1780=(ad1780-repmat(mean(ad1780(903:end,:)),981,1))./repmat(std(detrend(ad1780(903:end,:))),981,1);

load data1800.txt
ad1800=[zeros(800,102); data1800(:,2:end)];  
for ii=1:981
    iNaN=isnan(ad1800(ii,:));  %% if there are NaNs, replace with previous 
    if sum(iNaN)
        fiNaN=find(iNaN);
        ad1800(ii,fiNaN)=ad1800(ii-1,fiNaN);
    end
end


AD1800=(ad1800-repmat(mean(ad1800(903:end,:)),981,1))./repmat(std(detrend(ad1800(903:end,:))),981,1);

load data1820.txt
ad1820=[zeros(820,112); data1820(:,2:end)];

for ii=1:981
    iNaN=isnan(ad1820(ii,:));  %% if there are NaNs, replace with previous 
    if sum(iNaN)
        fiNaN=find(iNaN);
        ad1820(ii,fiNaN)=ad1820(ii-1,fiNaN);
    end
end

AD1820=(ad1820-repmat(mean(ad1820(903:end,:)),981,1))./repmat(std(detrend(ad1820(903:end,:))),981,1);

% Instrumental data singular value decomposition, A=U*S*Vt 

load pc01.txt
load pc02.txt
load pc03.txt
load pc04.txt
load pc05.txt
load pc06.txt
load pc07.txt
load pc08.txt
load pc09.txt
load pc10.txt
load pc11.txt
load pc12.txt
load pc13.txt
load pc14.txt
load pc15.txt
load pc16.txt


U=[pc01(:,2) pc02(:,2) pc03(:,2) pc04(:,2) pc05(:,2) pc06(:,2) pc07(:,2) pc08(:,2) pc09(:,2) pc10(:,2) pc11(:,2) pc12(:,2) pc13(:,2) pc14(:,2) pc15(:,2) pc16(:,2)];

load tpca_eigenvals.txt

S=diag(tpca_eigenvals(1:16,2)); % First 16 values

load eof01.txt
load eof02.txt
load eof03.txt
load eof04.txt
load eof05.txt
load eof06.txt
load eof07.txt
load eof08.txt
load eof09.txt
load eof10.txt
load eof11.txt
load eof12.txt
load eof13.txt
load eof14.txt
load eof15.txt
load eof16.txt

Vt=[eof01(:,2) eof02(:,2) eof03(:,2) eof04(:,2) eof05(:,2) eof06(:,2) eof07(:,2) eof08(:,2) eof09(:,2) eof10(:,2) eof11(:,2) eof12(:,2) eof13(:,2) eof14(:,2) eof15(:,2)  eof16(:,2)]';



load gridpoints.txt


i_nh=find(gridpoints(:,4)>=0); % indices for NH in i_nh

% stds from tpca_standard.txt
load tpca_standard.txt
sigmas=tpca_standard(i_nh,3);


% Reference temperature

fid1=fopen('nhmean.txt','r');
dds=[];
tline=fgetl(fid1);
while 1
            tline = fgetl(fid1);
            if ~ischar(tline), break, end
            dds=[dds;tline];
end

fclose(fid1);

ddn=str2num(dds);

RefT=ddn(503:581,3); 

% Sparse reference

fid1=fopen('nhem_sparse.dat','r');
dds=[];
for iz=1:34
    tline=fgetl(fid1);
end
while 1
            tline = fgetl(fid1);
            if ~ischar(tline), break, end
            dds=[dds;str2num(tline)];
end

fclose(fid1);

RefT_ver=dds(1:48,2); % Verification reference 1854-1901, sparse grid

SRefT=dds(49:127,2); % Sparse for calibration period


i_ver=855:902; % Verification indices for reconstruction
i_cal=903:981; % Calibration indices


% cosines of latitudes (Cos) , i_nh indices to NH
Cos=cos(gridpoints(i_nh,4)*pi/180);


t=(1000:1980)'; % reconstruction years

% Reconstruction stats, see http://www.climateaudit.org/?p=647
RE=[]; % calibration REs
RE_ver=[];
OUT=[]; % output
s2CAL=[];
R2=[];
R2_ver=[];
CE=[];
s2VER=[];


%%% PART III       %%%  
%%% Perform the calibration in 12 steps %%%


% Plot the original reconstruction
figure(1)
subplot(211)
plot(t,nhem_recon(:,2))
hold on
title('Original in blue, PC1 not fixed ')
xlabel('Year')
ylabel('Temperature') 

subplot(212)
plot(t,nhem_recon(:,2))
hold on
title('Original in blue, PC1 fixed ')
xlabel('Year')
ylabel('Temperature') 


for step_i=1:13 %% Run the calibration step by step, last for the pc1_fixed

    if step_i==1
    

    Pmc=AD1000; % calibration detrended variance =1, centered to cal mean
    
    
    
    Usel=[1]
    chkd=1:400;

    elseif step_i==2
    
        
        chkd=401:450;
    Pmc=AD1400;

   
    
    elseif step_i==3
     
    
    chkd=451:500;
    Usel=[1 2]  %% select U-vectors using http://holocene.meteo.psu.edu/shared/research/ONLINE-PREPRINTS/MultiProxy/stats-supp.html
    Pmc=AD1450;
    
   elseif step_i==4
       
       
    chkd=501:600;
    Usel=[1 2]
    Pmc=AD1500;
    
    
    
    elseif step_i==5

        
    
    chkd=601:700;
    Usel=[1 2 11 15]
    Pmc=AD1600;
    
    
    
    
    elseif step_i==6

    
    chkd=701:730;
    Usel=[1 2 5 11 15]
    Pmc=AD1700;
    
    
    
    elseif step_i==7

        chkd=731:750;
    Usel=[1 2 5 11 15]
    Pmc=AD1730;
    
    
    
    elseif step_i==8

    
    chkd=751:760;
    Usel=[1 2 3 5 6 8 11 15]
    Pmc=AD1750;
    
    
    elseif step_i==9

    
    chkd=761:780;
    Usel=[1 2 3 4 5 7 9 11 15]
    Pmc=AD1760;
    
    
    elseif step_i==10

    
    chkd=781:800;
    Usel=[1 2 3 4 5 7 9 11 14 15 16]
    Pmc=AD1780;
    
    
    elseif step_i==11

    
    chkd=801:820;
    Usel=[1 2 3 4 5 7 9 11 14 15 16]
    Pmc=AD1800;
    
    
    elseif step_i==12

        
    
    
    chkd=821:981;
    Usel=[1 2 3 4 5 7 9 11 14 15 16]
    
    % max RE with 11 Us
%Usel=[     1     2     3     4     5     7     8     9    11    12    13]
    
    
    Pmc=AD1820;
    
    else

     Pmc=AD1000b; % 
     Usel=[1]
     chkd=1:400;

        
    
    end % if step_i    

% http://www.climateaudit.org/?p=2445
zz=pinv(U(1:79,Usel))*Pmc(903:end,:);  % part one of CCE, temperature PCs as targets
zzz=Pmc(1:end,:)*pinv(zz); % part two of CCE




Uhscaled=zeros(length(zzz),length(Usel)); 

for zii=1:length(Usel) % variance matching, 

Uhscaled(:,zii)=zzz(:,zii)/std(zzz(903:end,zii))*std(U(1:79,Usel(zii)));

%Uhscaled(:,zii)=zzz(:,zii) % Don't touch the variance

end


temp_diag=diag(S);
temp_diag=temp_diag(Usel);
S_=diag(temp_diag); % Selected values of S diagonal


% http://www.climateaudit.org/?p=520
%
% http://www.climateaudit.org/?p=530

Recs=Uhscaled*S_*Vt(Usel,:); % back to temperatures 

RecsN=Recs(:,i_nh)*diag(sigmas)*ones(length(i_nh),1)/sum(Cos); % ?, but works

RecNH=RecsN; % 


    
V1=nhem_recon(chkd,2);
V2=RecNH(chkd,1);
    
    
    

figure(1)
if step_i<13
    
    OUT=[OUT; t(chkd) V2 V2]; % OUT= [time non-fixed fixed]

    subplot(211)
    plot(t(chkd),V2,'g')
    
    if step_i > 1
    subplot(212)
    plot(t(chkd),V2,'g')
    end
else
    subplot(212)
    plot(t(chkd),V2,'g')
    OUT(1:400,3)=V2;
end



Residual=RecNH(903:981)-RefT;

Residual_ver=RecNH(i_ver)-RefT_ver;




s2CAL=[s2CAL; t(chkd(1)) 2*std(Residual)]; % Residual zero mean

s2VER=[s2VER; 2*sqrt(Residual_ver'*Residual_ver/length(Residual_ver)) ];

RE=[RE; t(chkd(1)) 1-Residual'*Residual/(RefT'*RefT)];

ctemp=corrcoef(RecNH(903:981),RefT);

R2=[R2; t(chkd(1)) ctemp(2,1)^2];

RE_ver=[RE_ver; t(chkd(1)) 1-Residual_ver'*Residual_ver/(RefT_ver'*RefT_ver)]; % mean(RefT)=0

CE=[CE; 1-Residual_ver'*Residual_ver/((RefT_ver-mean(RefT_ver))'*(RefT_ver-mean(RefT_ver))) ];

ctemp=corrcoef(RecNH(i_ver),RefT_ver);

R2_ver=[R2_ver; t(chkd(1)) ctemp(2,1)^2 ];


end % for i=

format long g

disp('Calibration R2: ')
 R2
 
disp('Verification R2: ')
 R2_ver
 

 
% plot the AD 1000-1850 trends

XX=[ones(851,1) (1000:1850)'];
Tnf=pinv(XX)*OUT(1:851,2);
Tf=pinv(XX)*OUT(1:851,3);

% trends per century
Tnf_c=Tnf(2)*100;
Tf_c=Tf(2)*100;

Trendnf=XX*Tnf;
Trendf=XX*Tf;



figure(1)
subplot(211)
plot((1000:1850)',Trendnf,'r')
text(1500,0.25,['Trend per century: ', num2str(Tnf_c,3) ])
subplot(212)
plot((1000:1850)',Trendf,'r')
text(1500,0.25,['Trend per century: ', num2str(Tf_c,3) ])


% That's all