clear
n=40;
nt=25;
threshold = 0.2;
MC = 1000;
C = zeros(MC,1)+NaN;
for mc=1:MC
  x=randn(n,nt);
  y = randn(1,nt);
  %compute correlation between y and pool of predictors
  c = correlation(x,y(ones(n,1),:));
  %choose the ones whose correlation exceeds the threshold
  good = find(abs(c)>threshold);
  n2 = min([nt,length(good)]);
  x2 = x(good,:);

  %PCA of retained predictors
  [u,s,v] = svd(x2,'econ');
  x2 = v';

  %keep PCs whose corelation exceeds the threshold
  c2 = correlation(x2,y(ones(n2,1),:));
  good2 = find(abs(c2)>threshold);
  x3 = x2(good2,:);
    
  %cross-validated regression
  for ii=1:nt
    ind = setdiff([1:n],ii);
    b= regress(y',[x3;ones(1,nt)]');
    xp(ii) = [x3(:,ii);1]'*b;
  end
  C(mc)= corr2(y,xp);
end
return

plot(1982:2006,[y;xp]','-x')
title(['correlation =',num2str(corr2(y,xp),2)])
legend('obs.','CV forecast')
filename='screen_time_series';
export_fig('-pdf',filename)

hist(C)
title('frequency of correlation')
filename='screen_hist';
export_fig('-pdf',filename)
