% Exercise 13 - Working with increments

% Open the file Tsup.mat with the Lake Superior air temperatures.
% Quickly repeat the results of the previous exercise by running 
% this Matlab cell.
% You can execute this .m file one cell at a time by clicking on 
% the corresponding button in the toolbar.

t=Tsup(:,1);
T1=Tsup(:,2);
T2=Tsup(:,3);

plot(t,T1,t,T2)
%%
ws=24;   % size of averaging window is 1 day
rmean1=filter(ones(1,ws)/ws,1,T1);
Q1 = T1 - rmean1;

rmean2=filter(ones(1,ws)/ws,1,T2);
Q2 = T2 - rmean2;
% The new series Q1 and Q2 now have the annual cycle removed.

%%
plot(t,Q1,t,Q2)
%%
plot(Q1,Q2,'o')
%%
c=xcorr(Q1,Q2,'unbiased'); plot(c)
% If you zoom in on the graph you can see the daily cycle (24 hours).

%% 
% A nice trick in time series analysis is working with the increments,
% rather than with the series itself.
% Let's construct the series of hourly _increments_ in temperature. 
% That is, the i-th element of the increment series is 
%   P1(i) = T1(i+1)-T1(i)
% A simple programming loop will do this:
d=max(size(t));
for i=1:d-1
    P1(i)=T1(i+1)-T1(i);
    P2(i)=T2(i+1)-T2(i);
end
tinc(1:d-1)=t(1:d-1); % new time vector is 1 element shorter than t
%%
plot(tinc,P1,tinc,P2)
% Notice that the sinusoidal annual cycle is removed
% but annual variations are still noticable because
% hourly temperature _fluctuations_ are larger in some seasons.
% Compare this plot to the temperature record to find out
% which season that is.
%%
plot(P1,P2,'o')

% Notice that P1 and P2 are less correlated than Q1 and Q2 --
% chances in weather are less synchronous between the two locations
% than the weather itself!
%%
c=xcorr(P1,P2,'unbiased'); plot(c)

% Notice any cycles in the data.
% See if the daily cycle is still detectable.

%%