Code covered by the BSD License  

Highlights from
deseason

image thumbnail

deseason

by

 

26 Aug 2013 (Updated )

Remove the seasonal signal from a daily time series of data.

deseason(t,y)
function [yds,ymean,annualsignal] = deseason(t,y)
% [yds,ymean,annualsignal] = deseason(t,y) returns a time series yds of the
% same length as the multi-year daily time series y.  ymean is a vector the
% same length as t and y, containing only the seasonal signal.  The vector
% annualsignal is of length 366, the annual signal of daily values for day
% of year 1 to 366. Note that this will only work for daily records with t
% in datenum format.  
% 
% This function works by finding the mean value of the signal corresponding to 
% for each day of the year. That mean value is then subtracted from the raw signal 
% corresponding to that day of every year. For leap years, the annualsignal
% value of the 366th day of the year is calculated as a simple average of
% the annualsignal for Dec. 31 and Jan. 1. 
% 
% Please note that use of this function may be considered somewhat sloppy
% science. It's a clunky method of removing seasonal signals, but if your
% time series is long enough this deseason function is a quick and easy way
% to compare the magnitude of seasonal variations versus other variations
% in your signal. Or simply plotting annualsignal is a fair qualitative
% measure of the overall shape of a seasonal signal.  
%
% * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
% REQUIREMENTS: 
% 
% This function requires Anthony Kendall's date2doy function available on the 
% Mathworks FEX site: http://www.mathworks.com/matlabcentral/fileexchange/18316
% 
% * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
% EXAMPLE: Deseason some daily time series of temperature or something: 
% 
% % Declare some input data with seasonal signal in it: 
% t = (datenum('Jan 5, 2000'):datenum('May 15, 2016'))'; 
% temp = 25*cos(2*pi*date2doy(t)/365.25) + 30*rand(size(t))-15; 
% 
% % Deseason the data:
% [yds,ymean] = deseason(t,temp);
% 
% % Plot the raw and deseasoned data:
% figure
% set(gcf,'position',[100 100 600 250]); 
% plot(t,temp)
% datetick
% ylabel('temperature or something');
% axis tight
% box off
% hold on;
% plot(t,yds,'k')
% plot(t,ymean,'r')
% legend('raw daily data','deseasoned data','seasonal signal',...
%     'location','eastoutside'); 
% legend boxoff
% 
% * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
% Created by Chad A. Greene, August 2013. 
% version 2: now restores DC offset to yds, January 2014. 


% Input check:
if length(t)~=length(y)
    disp('Error: vectors t and y must be same length.')
    return
end

t = floor(t); % ensures we're using simple daily values 

tdoy = date2doy(t); 

annualsignal = NaN(366,1); 
yds = NaN(size(y)); % this'll be the deseasoned signal. 
ymean = NaN(size(y)); % This'll be a time series of seasonal signal the length of the full record. 

for d = 1:365
    annualsignal(d) = mean(y(tdoy==d)); 
    yds(tdoy==d) = y(tdoy==d)-annualsignal(d); 
    ymean(tdoy==d) = annualsignal(d); 
end

% Linearly interpolate a mean value for Dec 31 of leap years: 
annualsignal(366) = (annualsignal(365)+annualsignal(1))/2;
ymean(tdoy==366) = annualsignal(366); 
yds(tdoy==366) = y(tdoy==366)-annualsignal(366); 

% Restore the DC offset: 
yds = yds + mean(y); 

end

Contact us