function [lengths,growths,decays,event_dates,peak_strengths,mean_strengths] = eventLengths(time_series,date_vector,threshold,type)
% Takes a time series of anomalies, date vector, and threshold and finds
% all events and their lengths in the time series for either positive or
% negative events. Event lengths are given in units of time step
% of the time axis, assuming that time stepping is regular.
%
% USAGE: [lengths,growths,decays,event_dates,peak_strengths,mean_strengths] = ...
% eventLengths(time_series,date_vector,threshold,type)
%
% lengths, growths,
% decays, event_dates: vectors that contain one value per event.
% time_series: the data, a time series of anomalies
% threshold: threshold above the mean (0 for anomalies) for which
% an anomaly is potentially part of an event
% type: 1 for positive events, -1 for negative events
%
% An event is identified when anomalies of the same sign exceed a
% threshold with respect to the standard deviation of the anomaly type. E.g.,
% for positive events, take calculate the mean of all positive events
% and apply the threshold for identification to this mean,
% instead of to the standard deviation of the entire anomaly
% time series. (--> This is what eventLengths() does.)
%
% In comparison to eventLengths(), relativeEventLengths() sets
% the threshold for events differently and allows for
% systematic differences between positive and negative
% events. This is especially important when tackling
% issues of ENSO asymmetry.
% ------------- %
% QUALITY CHECK %
% ------------- %
% date vector and data of same size?
n_ts = length(time_series);
n_time = size(date_vector,1);
if ( n_ts ~= n_time )
error('Dimension missmatch between time series and date vector.');
end
% ------------ %
% PREPARE DATA %
% ------------ %
% flip data upside down for negative events
anomalies = type.*time_series;
% mean of anomalies of interest, i.e. of either positive or
% negative anomalies; for type = -1, the time series has
% been flipped already
% identify all desired anomalies
is_anomaly = ( anomalies > 0 );
relative_base = std(anomalies(is_anomaly));
% for verification purposes
complementary_base = std(anomalies(~is_anomaly));
% move offset according to threshold
offset = anomalies-(threshold*relative_base);
% date numbers for easier handling
date_numbers = datenum(date_vector);
% --------------- %
% IDENTIFY EVENTS %
% --------------- %
% set the time series to 1 for all positive events
base_line = sign(offset);
base_line(base_line<0) = 0;
% breaks between events
event_boundaries = diff(base_line);
is_boundary = ( event_boundaries ~= 0 );
% find indivudal events
event_starts = find(event_boundaries > 0);
event_stops = find(event_boundaries < 0);
% if there is either no event or no event_stops:
% return NaN values
if ( numel(event_starts) == 0 | numel(event_stops) == 0 )
lengths = NaN;
growths = NaN;
decays = NaN;
event_dates = NaN;
peak_strengths = NaN;
mean_strengths = NaN;
return;
end
% ignore partial events at the beginning nd end of the time series
% leading events: ignore the first partial event
if event_stops(1) < event_starts(1)
event_stops = event_stops(2:end);
end
% trailing event: ignore the last (started) partial event
if event_starts(end) > event_stops(end)
event_starts = event_starts(1:end-1);
end
% add +1 to all event starts to make sure that they only
% start when the anomaly value is already higher than the threshold
event_starts = event_starts+1;
% find events that last only for 1 month and ignore them
micro_events = find( (event_stops-event_starts) < 2 );
event_starts(micro_events) = [];
event_stops(micro_events) = [];
% a check: there must be as many starting events as ending events
if ( length(event_starts) ~= length(event_stops) )
error('Missed leading or trailing events.');
end
% -------------------------------- %
% IDENTIFY PEAKS AND COUNT LENGTHS %
% -------------------------------- %
% number of events
n_events = length(event_starts);
% safe data
lengths = zeros(1,n_events);
growths = zeros(1,n_events);
decays = zeros(1,n_events);
event_dates = zeros(1,n_events); % as datenums
peak_strengths = zeros(1,n_events);
mean_strengths = zeros(1,n_events);
% loop over all events
for c_event = 1:n_events
% IDENTIFY EVENTS AND DIAGNOSE THEIR LENGTHS
% ++++++++++++++++++++++++++++++++++++++++++
% event cropping
c_start = event_starts(c_event);
c_stop = event_stops(c_event);
event_data = anomalies(c_start:c_stop);
% total event length
event_length = length(c_start:c_stop);
lengths(c_event) = event_length;
% growing and decaying phases (this is a bit tricky,
% tread cautiously here)
peak_position = find(event_data == max(event_data));
% if more than one peak exists: good chance of corrupted data
% ignore
if ( length(peak_position) > 1 )
lengths = NaN;
growths = NaN;
decays = NaN;
event_dates = NaN;
peak_strengths = NaN;
mean_strengths = NaN;
return;
end
% absolute position for the date vector
abs_peak = c_start+peak_position-1;
% save date of peak
event_dates(c_event) = date_numbers(abs_peak);
% growth and decay lengths in months
% count peak-position twice, since growth ends there
% and decay starts
growths(c_event) = length(1:peak_position);
decays(c_event) = length(peak_position:event_length);
% event strengths: peak strength and mean strength
peak_strengths(c_event) = type*anomalies(abs_peak); % to restore the correct sign
mean_strengths(c_event) = type*mean(event_data);
end
% convert event peak date numbers into a more usable date vector
event_dates = datevec(event_dates);