% autocorrelation:
rho = 0.95;

n=2000;

rr = randn(n,1000);
for k = 1:n
  disp(k);
  fflush(stdout);
  % spin-up:
  x(1,k)=0;
  for j=1:100
    x(1,k)=rho*x(1,k)+(1-rho)*rr(k,j);
  end
  % Gauss-Markov:
  for i = 2:1000
    x(i,k) = rho*x(i-1,k)+(1-rho)*rr(k,i);
  end
end

% Filter out trends, trend function
for i=1:100
  tf(i) = (i-100)/50+1;
end

for k = 1:n
  c = cov(x(901:1000,k),tf');
  v1= cov(x(901:1000,k));
  v2=cov(tf);
  % correlation:
  c = c/sqrt(v1*v2);
  % 'screening':
  if (c < 0.6)
    for i=1:1000
      x(i,k) = NaN;
    end
  else
    j=k;
  end
end

printf('valid members: %d\n\n', n-sum(isnan(x(1,:))));

for i = 1:1000
  a=x(i,:);
  f(i)=quantile(a',0.25);
  m(i)=quantile(a',0.5);
  t(i)=quantile(a',0.75);
end

plot(1001:2000,x(1:1000,j),'-b',1001:5:2000,m(1:5:1000),'-r',1001:5:2000,f(1:5:1000),':m',1001:5:2000,t(1:5:1000),':m', 1901:2000,tf,'-g',[1001 2000],[0 0],'-k');
legend('Example', 'median', '25% percentile', '75% percentile', 'trend function', '', 'location', 'northwest');

print('hockey.eps', '-depsc');
print('hockey.png', '-dpng');