clearvars -except fac fac2
nt = 10000;
n=12
q0 = zeros(n,n);q0=eye(n);
q = fac2*q0 +  (1-fac2)*ones(n);
a = [[zeros(n-1,1),eye(n-1)];zeros(1,n)];

a = zeros(n);
a1 = diag(diag(eye(n)),1);
a2 = diag(diag(eye(n)),2);
a = a1(1:n,1:n);


sh = a;
sh(n,1)=0;
a = sh;

a = a + fac*eye(n);
a = a/max([svd(a)])*0.9;
%a(logical(eye(n)))=0.2;
%v = randn(n);
%[q,r] = qr(v);
%a = q*diag(rand(1,n))*q';

x = zeros(n,nt);
x(:,1)=0;
for ii=1:(nt-1)
  x(:,ii+1) = a*x(:,ii) + q*randn(n,1);
end

x(:,1)=x(:,nt);

for ii=1:(nt-1)
  x(:,ii+1) = a*x(:,ii) + q*randn(n,1);
end

k=12;
f = NaN(k,n,nt+k);
for ii=1:k
  wrk = a^ii;
  f(ii,:,(ii+1):(nt+ii)) = wrk*x;
end

subplot(3,1,1)
my_pcolor([1:50],[0:k],[x(1,1:50);squeeze(f(:,1,1:50))]);shading('faceted')
cm = 0.9*max(abs(f(:)));
caxis([-cm,cm])
xlabel('T')
ylabel('L')
set(gca,'ytick',[0:2:12],'yticklabel',{'Obs.',int2str([2:2:12]')})
oo = x(ones(k,1),(k+1):nt);
ff = squeeze(f(:,1,(k+1):nt));

nlags = 9;
for ii=1:k
    lc(ii,:) = xcorr(oo(ii,:),ff(ii,:),nlags,'coeff');
end

subplot(3,1,2)
plot([-nlags:nlags],lc(1:2:9,:))
caxis([-0.3,1])
title('lag correlation')
grid
legend('1','3','5','7','9')
set(gca,'ylim',[-0.05,1.05])
hold on 
plot([0,0,],[0,1],'linewidth',2,'color','k')
[yy,idum]=max(lc(1:2:9,:)');
xx = [-nlags:nlags];
plot(xx(idum),yy,'ko')
plot(xx(idum),yy,'k+')
hold off

C = q;
for ii=1:10
C = C + a*C*a';
end

qq = eye(n);
qq = qq(:,1);
qq = qq/sqrt(qq'*C*qq);
d1=0;
for ll = -nlags:nlags
d1=d1+1;
  d2=0;
  for kk =1:2:9
    d2=d2+1;
    xl = x(:,(1+nlags+9-ll):end);
    xk = x(:,(1+nlags+9-kk):end);
    m = min([length(xk),length(xl)]);
    cc(d1,d2) = corr2(qq'*a^kk*xk(:,1:m),qq'*xl(:,1:m));
    cc2(d1,d2) = (qq'*a^kk*C*(a')^abs(kk-ll)*qq)/sqrt(qq'*a^kk*C*(a')^kk*qq);
  end
end
subplot(3,1,3)
plot([nlags:-1:-nlags],cc)
plot([nlags:-1:-nlags],cc2)
return

nlags = 9;
for ii=1:n
    lc_another(ii,:) = xcorr(x(ii,:),x(1,:),nlags,'coeff');
end

pcolor([n:-1:1],[-nlags:1:nlags],lc_another')
xlabel('"basin"')
ylabel('lag')
title('lag correlation across basin')
caxis([0,1])
