% Model:
% Temp = a*CO2 + b*Solar+c*Volcanic+d+n
% n is i.i.d Gaussian, std ~ 0.1
% 
%

load fig7-co2.txt;
load fig7-corrs.txt;
load fig7-nh.txt;
load fig7-solar.txt;
load fig7-volcanic.txt;

window=200;

first=1;
last=window;

OUT=[];

zz=1000;
yr=fig7_nh(window/2,1);
for i=1:zz
    
    t_co2=fig7_co2(first:last,2);
    std_co2=std(t_co2);
    
    t_solar=fig7_solar(first:last,2);
    std_solar=std(t_solar);
    
    t_volcanic=fig7_volcanic(first:last,2);
    std_volcanic=std(t_volcanic);
    
    t_nh=fig7_nh(first:last,2);
    std_nh=std(t_nh);
    
    % LS solution 
    A=[t_co2 t_solar t_volcanic ones(window,1)];    
    out=inv(A'*A)*A'*t_nh;
    
    out=(out/std_nh).*[std_co2; std_solar; std_volcanic;1];
    
    OUT=[OUT; yr out']; % yes, yes, this is the slow way.. 
    
    if last+1 <= length(fig7_nh)
        first=first+1;
        last=last+1;
        yr=yr+1;
    else
        break
    end
end

 close all
    plot(OUT(:,1),OUT(:,2:4))
    legend('CO2','Solar','Volcanic')
    hold
    plot(fig7_corrs(:,1),fig7_corrs(:,2:4),'-.')
    
    % use all data
    Aa=[fig7_co2(:,2) fig7_solar(:,2) fig7_volcanic(:,2) ones(length(fig7_co2),1)];
    outa=inv(Aa'*Aa)*Aa'*fig7_nh(:,2);
	 repl=Aa*outa;
    
    
    % simulated temps, assume repl is correct
    randn('state',sum(100*clock))
    n=randn(length(fig7_co2),1)/10;
    s_fig7_nh=repl+n;
    repls=Aa*inv(Aa'*Aa)*Aa'*s_fig7_nh;
    
    figure
    subplot(211)
    plot(repl)
    hold
    plot(fig7_nh(:,2),'r.')
    legend('Adjusted','Original')
    subplot(212)
    plot(repls)
    hold
    plot(s_fig7_nh,'r.')
    legend('Adjusted_S','Original_S')
    
    figure
    residual=repl-fig7_nh(:,2);
    residuals=repls-s_fig7_nh;
    subplot(211)
    plot(residual)
    legend('Residual')
    std(residual)
    subplot(212)
    plot(residuals)
    legend('Residual_S')
    std(residuals)
    
    
    figure
    fft_res=fft(residual);
    fft_ress=fft(residuals);
    subplot(211)
    plot(fft_res.'.*fft_res')
    legend('PSD (not scaled)')
    subplot(212)
	 plot(fft_ress.'.*fft_ress')
    legend('PSD_S (not scaled)')
    

    Source: geocities.com/uc_edit