% MBH98 Figure 7 emulation, 
% ref http://climateaudit.org/2006/05/31/more-on-mbh98-figure-7/ 
% UC 2011

load fig7-co2.txt;  % http://www.nature.com/nature/journal/v430/n6995/extref/FigureData/fig7-co2.txt
load fig7-corrs.txt; % http://www.nature.com/nature/journal/v430/n6995/extref/FigureData/fig7-corrs.txt
load fig7-nh.txt; % http://www.nature.com/nature/journal/v430/n6995/extref/FigureData/fig7-nh.txt
%  ( comment the first three lines of fig7-nh.txt )

load fig7-solar.txt; % http://www.nature.com/nature/journal/v430/n6995/extref/FigureData/fig7-solar.txt
load fig7-volcanic.txt; % http://www.nature.com/nature/journal/v430/n6995/extref/FigureData/fig7-volcanic.txt

window=201; % this gives the best agreement with fig7_corrs

%window=100 % MBH98: Nonetheless, all of the important conclusions drawn
            %         below are robust to choosing other reasonable (for example, 100-
            %         year) window widths.

first=1;
last=window;

n_i=length(fig7_co2)-window+1;

OUT=zeros(n_i,7);   
zz=n_i; %

yr=fig7_nh(round(window/2),1);

% AR(1) covariance matrix
rho=0.77; % MBH98
rho_v=rho.^(0:window-1);
S=toeplitz(rho_v);

% Averaging matrix
M1=ones(window)./window;
M=eye(window)-M1;


% 90, 95, 99 for normal (norminv(0.9,0,1)..)
L=[1.2816 1.6449 2.3263]';


for i=1:zz
    
    t_co2=(fig7_co2(first:last,2));
    m_co2=mean(t_co2); 
    std_co2=std(t_co2);
    t_co2=(t_co2-m_co2)/std_co2; % average & scale co2
    
    t_solar=(fig7_solar(first:last,2));
    m_solar=mean(t_solar); 
    std_solar=std(t_solar);
    t_solar=(t_solar-m_solar)/std_solar; % average & scale solar   
    
    
    
    t_volcanic=(fig7_volcanic(first:last,2));
    m_volcanic=mean(t_volcanic); 
    std_volcanic=std(t_volcanic);
    t_volcanic=(t_volcanic-m_volcanic)/std_volcanic; % average & scale volcanic     
    
    t_nh=standardize(fig7_nh(first:last,2));
    m_nh=mean(t_nh); 
    std_nh=std(t_nh);
    t_nh=(t_nh-m_nh)/std_nh; % average & scale nh reconstruction     

    
    
    % LS solution; nh,co2,solar,volcanic now standardized 
    A=[t_co2 t_solar t_volcanic ];  
    
    pA=pinv(A);
    out=pA*t_nh;
    
    
    % Check what would Gaussian AR(1) 'nh reconstruction' give 
    % (**Note: unity variance process, not
    % standardized, will do to prove the point**)
    % MBH98: "The associated
    %         confidence limits are approximately constant between sliding 200-
    %         year windows."
             
             
    S_A=pA*M*S*M'*pA';
    
    L_CSV=sqrt(diag(S_A));
    
    OUT(i,:)=[ yr out' L_CSV']; % store 
    
        first=first+1; % move the window one year
        last=last+1;
        yr=yr+1;
end

close all

% plot and compare with Nature Figure Data
plot(OUT(:,1),OUT(:,2:4),'o')
legend('CO2','Solar','Volcanic')
hold

if window==201
plot(fig7_corrs(:,1),fig7_corrs(:,2:4),'x')
end

plot(OUT(:,1),OUT(:,5)*L(3))
plot(OUT(:,1),OUT(:,5)*L(2),'--')
plot(OUT(:,1),OUT(:,5)*L(1),':')


plot(OUT(:,1),OUT(:,6)*L(3),'g')
plot(OUT(:,1),OUT(:,6)*L(2),'g--')
plot(OUT(:,1),OUT(:,6)*L(1),'g:')


plot(OUT(:,1),-OUT(:,7)*L(3),'r')
plot(OUT(:,1),-OUT(:,7)*L(2),'r--')
plot(OUT(:,1),-OUT(:,7)*L(1),'r:')


plot(fig7_corrs(:,1),fig7_corrs(:,5:7),'k--')

if window==201
figure
plot(round(OUT(:,2:4)*1000)/1000-fig7_corrs(:,2:4)) % round & check 
end