% Exercise 14 - Significance of autocorrelation. Delays. Detection of cycles.

% Open the file GreenAntarc.mat.
% The file contains climate data from ice cores in Greenland and Antactica
% The first column is age in thousands of years (kyr) before present.
% The second and third colums are the concentrations of oxygen-18
% isotope in Greenland and Antarctica, respectively.
% The 18-O concentration is a proxy for the regional temperature,
% so higher concentrations mean warmer climate.

% Here, we will look at correlations and evaluate shifts in the data
% The goals is to answer how climate changes in Greenland
% are related to those in Antarctica.

% Copy data to separate variables. Notice that the time
% is sampled at 0.1 kyr intervals.

t=GreenAntarc(:,1);
TG=GreenAntarc(:,2);
TA=GreenAntarc(:,3);
plot(t,TG,t,TA)

% Visually inspect the series. Are there any apparent cycles? Shifts?
% Do the series appear stationary in time?
%%
% Calculate the running means for both series
ws=20;   % <--- Experiment with size of the averaging window (in 0.1 kyr).
rmeanG=filter(ones(1,ws)/ws,1,TG);
rmeanA=filter(ones(1,ws)/ws,1,TA);
plot(t,rmeanG,t,rmeanA)
%%
% This cell is necessary to remove some artifacts of averaging
for i=1:100
    if rmeanG(i)>-34
        rmeanG(i)=TG(i);
    end
    if rmeanA(i)>-34
        rmeanA(i)=TA(i);
    end
end
plot(t,rmeanG,t,rmeanA)
%%
% Subtract the running mean from the original series
sG=TG-rmeanG;
sA=TA-rmeanA;
plot(t,sG,t,sA)
% The new series sG and sA now have only high-frequency "noise"
%%
% Now compute the correlations and autocorrelations
% Modify this code to suit your task.
% Compute x-correlations, autocorrelations, correlations for running means.
% Experiment with different sizes of averaging window.

c=xcorr(sA,'coeff');plot(c)   % <-- autocorrelation
%c=xcorr(sG,sA,'coeff');plot(c) % <-- x-correlation
hold on
mc=0.5;
d=801*2-1;  % The series contains 801 points

plot([d/2 d/2],[-mc mc],'r')
hold off
%%
% You must have noticed that the series appear shifted
% relative to each other.
% Plot correlation coefficients for different shift values
for j=13:18    % < -- modify the range of shifts (in units of 0.1 ky)
sh=j;
G3=NaN(801-sh,1);  
A3=NaN(801-sh,1);
t3=NaN(801-sh,1);
for i=1:801-sh
    t3(i)=t(i);
    
    G3(i)=TG(i);         
    A3(i)=TA(i+sh);   % this shifts the original Antarctic series
    %G3(i)=sG(i);
    %A3(i)=sA(i+sh);   % this shift the "noise-only" Antarctic series
end
sh          % print shift value to screen
r=corrcoef(G3,A3);
r(2,1)      % print correlation coefficient
end
% Look at the Matlab output. Write down the maximum correlation
% coefficient value. Re-run this cell choosing the value of j that 
% corresponds to the shift that produces maximum correlation.

%%
plot(t3,G3,t3,A3)
%%
plot(G3,A3,'o')
%%
% plot the correlations for shifted series 
% Notice that the value of shift (variable 'sh') is as it was at the 
% end of the cycle above. Modify the code to plot the disired shift.

% Experiment by plotting different correlations and autocorrelations.
% Can you detect any cycles? Write down their periods!

c=xcov(G3,'coeff');plot(c)  
% c=xcov(G3,'coeff');plot(c)

% function 'xcov' works just like 'xcorr', only it first subtracts 
% the mean of the series.

hold on
mc=0.5;
d=(801-sh)*2-1;

plot([d/2 d/2],[-mc mc],'r')
hold off
%%
%Now let's determine if any of the correlations are significant.
% There is a test for this, similar to the t-test that we used before.
% Using the correlation coefficient 'r' at a given lag 'tau', we can 
% calculate the probability that 'r' is significantly different from zero.

r=r(2,1);    % <-- this uses corrcoef from the last tried shift
N=801-sh;
tau=sh;      % < --- assign the value of lag here
Z=r*sqrt(N-sh+3)
% The quantity Z is expected to be distributed Normally.
% A large Z means that the correlation is significant.
% The probability of obtaining a larger value of Z from uncorrelated
% data can be obtained from the CDF of a Normal distribution.