clear all; close all;

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

% the following two lines are for calculating
% with respect to the true construction
% (uncomment if needed)
%load nhmean.txt;
%nhmean(end-15:end,2)=fig7_nh(end-15:end,2);

n=386;
% window windth, change this to get interesting graphs!
W=200;

cc=0;
for k=0:n-W,
   window_nh=fig7_nh(k+1:k+W,2);
   % uncomment the line below if you want
   % with respect to true construction
   %window_nh=nhmean(k+211:210+k+W,2);
   window_co2=fig7_co2(k+1:k+W,2);
   window_solar=fig7_solar(k+1:k+W,2);
   window_volcanic=fig7_volcanic(k+1:k+W,2);

   % The following is "a direct way of
   % calculating"
   % the partial correlation
   % coefficients for the situation in hand
   % (see the appendix
   % of the link I provided over
   % ClimateAudit). You should have
   % this function directly in R.
   % (it is also in the statistical
   % toolbox in Matlab, which I don't
   % have access right now).
   % The ,1 thing in cov means
   % that Matlab normalizes with
   % respect to the length, not length-1.
   x=[window_nh window_co2]';
   y=[window_solar window_volcanic]';
   a=cov([x;y]',1);
   b=a(1:2,1:2)-a(1:2,3:4)*inv(a(3:4,3:4))*a(3:4,1:2);
   pcc(k+1,1)=b(1,2)/sqrt(b(1,1)*b(2,2));

   x=[window_nh window_solar]';
   y=[window_co2 window_volcanic]';
   a=cov([x;y]',1);
   b=a(1:2,1:2)-a(1:2,3:4)*inv(a(3:4,3:4))*a(3:4,1:2);
   pcc(k+1,2)=b(1,2)/sqrt(b(1,1)*b(2,2));

   x=[window_nh window_volcanic]';
   y=[window_solar window_co2]';
   a=cov([x;y]',1);
   b=a(1:2,1:2)-a(1:2,3:4)*inv(a(3:4,3:4))*a(3:4,1:2);
   pcc(k+1,3)=b(1,2)/sqrt(b(1,1)*b(2,2));
end
T=[(1609+W/2):(1995-W/2)];
plot(T,pcc)
legend('Co2','Solar','Volcanic',0);
hold on;
% The originals
plot(fig7_corrs(:,1),fig7_corrs(:,2:4),'-.');