% years 1861-1960

%load solarprox % solar data (solar), proxy data for (prox), nh recon (nh) 1861-1960


Rs=zeros(415,1);
Rn=Rs;

t=(1861:1960)';


for i=2:416

    r=corrcoef(solar(:,2),prox(:,i));
    Rs(i-1)=r(2,1);
    
    
    r=corrcoef(nh(:,2),prox(:,i));
    Rn(i-1)=r(2,1);
    
    
end

close all


proxc=prox;
proxc(:,254)=prox(:,254)*0+randn(100,1); % too redundant

DM=proxc(:,2:end);

bhat=(DM'*DM)\DM'*nh(:,2);


[notused,ST]=sort(Rs);
ST=flipud(ST)+1;
DMs=[proxc(:,ST(1:20)) ones(100,1)];

bhats=(DMs'*DMs)\DMs'*solar(:,2);

figure
plot(t,DMs*bhats)
hold on
plot(t,solar(:,2),'r')
legend('Solar reconstruction','Solar Original')



[notused,STn]=sort(Rn);
STn=flipud(STn)+1;
DMsn=[proxc(:,STn(1:20)) ones(100,1)];

bhatsn=(DMsn'*DMsn)\DMsn'*nh(:,2);


figure
plot(t,DMsn*bhatsn)
hold on
plot(t,nh(:,2),'r')
legend('NH reconstruction','NH Original')

    Source: geocities.com/uc_edit/mbh98

               ( geocities.com/uc_edit)