%%%% MBH99 residuals HOWTO
%%% ref http://www.climateaudit.org/?p=647
%% this is just a memo for myself, but can be useful for others as well
% uc_edit at yahoo.com


%%% get the data
% original reconstruction
load mbh99r.txt %in this subdirectory (source:
% http://www.clim-past-discuss.net/2/1001/2006/cpd-2-1001-2006-supplement.zip

t=(1000:1974)';

% proxies
load  mbh99p.txt %in this subdirectory (source:
% http://www.climateaudit.org/data/MBH99/proxy.txt
% note that pc1 is fixed using mbh99fix.dat, data years 1000-1974 only

% sparse temperature 
load sparse.txt % source http://www.ncdc.noaa.gov/paleo/ei/ei_data/nhem-sparse.dat
% years 1902-1975

% reverse engineer:
P=[mbh99p(1:400,:) ones(400,1)];
Ts=mbh99r(1:400,2);
betas=inv(P'*P)*P'*Ts;
mbh99rN=[mbh99p ones(975,1)]*betas; % how do we know that this works?

% we can never be sure, but here is a hint:
load MBH99_AD1000_res1.txt % scanned by Jean S from MBH99 submission
res=MBH99_AD1000_res1;
resN=mbh99rN(903:end)-sparse;

close all
plot(t,mbh99r(1:975,2))
hold on
plot(t,mbh99rN,'r')
legend('Original Full Rec','reverse eng. AD1000')

figure
plot(res)
hold on
plot(resN,'r')
legend('scanned','reverse eng., sparse ref')
corrcoef(res(1:73),resN)

    Source: geocities.com/uc_edit/MBH99

               ( geocities.com/uc_edit)