% cacsstsignificance.m 
%
% Demonstrate estimated the signficance levels for a set of EOFs

% The EOF analysis of CAC (Reynolds) SST data have to have been done by 
% running script cacsst.m

if exist('D')~=1
   disp('Running script cacsst to do the EOF analysis')
   pauseit = 0;
   cacsst
end
pauseit = 1;

% Plot fractional variance of first K modes
K = 15;
figure
plot(diag(S*S')/trace(S*S'),'--d')
set(gca,'xlim',[1 K])
legend('fraction of variance for each EOF mode')
xlabel('EOF mode number')
grid on

pause_it

% ----------------------------------------------------------------------
% Calculate Monte Carlo estimate of 95% significance level for a random
% noise data set of the same size as D
% ----------------------------------------------------------------------

[M N] = size(D);
nsamples = 10; % this is a bit small, but otherwise this step takes too long in class
f95 = eofsig_montecarlo(M,N,nsamples);

frac = diag(S*S')/trace(S*S');
plot(1:K,frac(1:K),'--d',1:K,f95(1:K),'r-d')
grid on
legend('fraction of variance for each EOF mode','95% significance level')
xlabel('EOF mode number')

pause_it

% For EOFs of data sets that have significant autocorrelation (say in
% the time dimension of the data matirx), it is recommended that N be
% reduced to Nstar = N/autocorrelation_scale to reflect the number of
% independent observations of the data in each row of the data matrix. 
% See Joliffe (2002) 'Principal Component Analysis,' Springer Series 
% in Statistics, chapter 6 and references therein. 

% We can see the effect by recomputing the EOFs for a decimated data set.

% Re-compute the EOFs using only every n'th value in the time series
figure
set(gcf,'pos',[18    39   521   646])
clf

n = 3; % decimate the time series by this factor

[U,S,V] = svd(D,0);
A = S*V';
frac = diag(S*S')/trace(S*S');
lab{1} = ['fraction \lambda_1 = ' num2str(frac(1),3)];
lab{2} = ['fraction \lambda_2 = ' num2str(frac(2),3)];
lab{3} = ['fraction \lambda_3 = ' num2str(frac(3),3)];
subplot(311)
pcolorjw(lon,lat,reshape(U(:,1),[nlats nlons]))
cax= caxis;
colorbar
text(145,15,lab,'color','w')
title('Original data set')

Dn = D(:,1:n:end);
[Un,Sn,Vn] = svd(Dn,0);
An = Sn*Vn';
fracn = diag(Sn*Sn')/trace(Sn*Sn');
lab{1} = ['fraction \lambda_1 = ' num2str(fracn(1),3)];
lab{2} = ['fraction \lambda_2 = ' num2str(fracn(2),3)];
lab{3} = ['fraction \lambda_3 = ' num2str(fracn(3),3)];
subplot(312)
pcolorjw(lon,lat,reshape(Un(:,1),[nlats nlons]))
caxis(cax)
colorbar
text(145,15,lab,'color','w')
title(['Data reduced in time by a factor of ' int2str(n)])

subplot(313)
pcolorjw(lon,lat,reshape(U(:,1)-Un(:,1),[nlats nlons]))
colorbar
pause
plot(year,A(1,:),year(1:n:end),An(1,:))
legend('a_1(t) original',['a_1(t) for every ' int2str(n) '^{th} data point']);
legend boxoff

pause_it

% The EOFs look the same because the monthly data resolve the variability
% very well. 

% But recompute the Monte Carlo significance estimate

[M N] = size(Dn);
nsamples = 10; % this is a bit small, but otherwise this step takes to long in class
f95_n = eofsig_montecarlo(M,N,nsamples);

figure
plot(1:K,frac(1:K),'--bd',1:K,f95(1:K),'b-d',...
   1:K,fracn(1:K),'--rd',1:K,f95_n(1:K),'r-d')
grid on
legend('fraction of variance','95% significance level',...
   ['fraction of variance for data reduced by ' int2str(n)],...
   '95% significance level')
xlabel('EOF mode number')

pause_it

% So how do we objectively estimate the effective degrees of freedom of the
% time sampling if the data? The answer lies in calculating the
% autocorrelation time scale, and assuming that the number of independent
% time samples is the duration of the time series divided by the
% autocorrelation time scale.

% One way of estimating the autocorrelation time scale is to form an
% ensemble autocorrelation function by averaging the auto-correlation 
% sequences estimated for each time series, then integrating this out to
% the first zero crossing i.e. the time lag at which, on average, a time 
% series beocmes uncorrelated with itself. 

% ----------------------------------------------------------------------
% Use my function degreesfreedom.m to estimate the effective d.o.f. 
% ----------------------------------------------------------------------
[dof,Tscale,xc,lags] = degreesfreedom(D);

disp(' ')
disp(['Estimated D.O.F. are ' int2str(dof)])
disp(' ')
disp(['Autocorrelation Tscale is ' num2str(Tscale) ' sampling ',...
      'intervals of the data'])

% Recompute significance limits using DOF = N/Tscale instead of N.
nsamples = 10; % this is a bit small, but otherwise this step takes to long in class
f95_dof = eofsig_montecarlo(M,dof,nsamples);

plot(1:K,frac(1:K),'-bd',1:K,f95(1:K),'b-s',...
   1:K,fracn(1:K),'-md',1:K,f95_dof(1:K),'m-s');
grid on
legend('fraction of variance','95% significance level',...
   'fraction of variance for data reduced by Tscale','95% significance level')
xlabel('EOF mode number')