Code covered by the BSD License  

Highlights from
Measures of Analysis of Time Series toolkit (MATS)

image thumbnail

Measures of Analysis of Time Series toolkit (MATS)


Dimitris Kugiumtzis


10 May 2010 (Updated )

MATS computes many measures of scalar time series analysis on many time series in one go.

function nuT = CorrelationDimension(xV,tauV,mV,theiler,sV,resol)
% nuT = CorrelationDimension(xV,tauV,mV,theiler,sV,resol)
% CORRELATIONDIMENSION computes the correlation dimension for a given time
% series 'xV', for a range of delays in 'tauV', a range of embedding 
% dimensions in 'mV' and for a range of upper/lower ratio of scaling window
% in % 'sV' (s=r2/r1 where r2-r1 is the length of the scaling window).
% The parameter 'theiler' excludes temporally close points (smaller than 
% 'theiler') from the inter-distance computations. The parameter 'resol'
% determines the number of radii for which the correlation sum is
% computed. 
% First, the correlation sum C(r) and the local slopes log(C(r))/log(r) are 
% computed for a range of distances r given by 'resol'. Then the correlation 
% dimension 'nu' is estimated by searching for the local slope in radii
% intervals [r1,r2] (determined by 's') with the smallest standard 
% deviation (best scaling). The mean local slope in this interval is the  
% estimate of 'nu'. 
% - xV      : Vector of the scalar time series ('xV' is then standardized
%             in [0,1]). 
% - tauV    : A vector of the delay times.
% - mV      : A vector of the embedding dimension.
% - theiler : the Theiler window to exclude time correlated points in the
%             search for neighboring points. Default=0.
% - sV      : A vector of values of upper/lower ratio of scaling window
%             (e=r2/r1 where r2-r1 is the length of the scaling window).
% - resol   : The number of radius for which the correlation sum is computed.
%             Note that this parameters is supposed to be larger than 10.
% - nuT     : A matrix of size 'ntau' x 'nm' x 'ne', where 'ntau' is the 
%             number of given delays, 'nm' is the number of given embedding 
%             dimensions and 'ne' is the number of scaling ratio of radii. 
%             The components of the matrix are the correlation dimension
%             values. 
%     <CorrelationDimension.m>, v 1.0 2010/02/11 22:09:14  Kugiumtzis & Tsimpiris
%     This is part of the MATS-Toolkit

% Copyright (C) 2010 by Dimitris Kugiumtzis and Alkiviadis Tsimpiris 
%                       <>

% Version: 1.0

%     This program is free software; you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation; either version 3 of the License, or
%     any later version.
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     GNU General Public License for more details.
%     You should have received a copy of the GNU General Public License
%     along with this program. If not, see>.

% Reference : D. Kugiumtzis and A. Tsimpiris, "Measures of Analysis of Time Series (MATS): 
% 	          A Matlab  Toolkit for Computation of Multiple Measures on Time Series Data Bases",
%             Journal of Statistical Software, in press, 2010

% Link      :
smoothorder = 2; % Order of the smooth derivative for the correlation sum
logdmin=-3; % the lowest log10(distance)
% Rescale to [0,1]
xmin = min(xV);
xmax = max(xV);
xV = (xV - xmin) / (xmax-xmin);
n = length(xV);
if isempty(theiler), theiler=0; end
nm = length(mV);
mV = sort(mV);
mmax = mV(end);
mmin = mV(1);
ntau = length(tauV);
ns = length(sV);
% Any computed distance d will be placed in a bin of 'rV'
rV = logspace(logdmin,0,resol);
rV = [-Inf rV]';
nuT = NaN*ones(ntau,nm,ns);
% To be used in allocating the components of temporal correlated distances
% in the distance matrix below (it is given here to avoid repetitions).

if resol-1-2*smoothorder < 1
if theiler>0
    isumV = cumsum([0:n-theiler-1]);
for itau = 1:ntau
    tau = tauV(itau);
    densityM = zeros(resol,mmax-mmin+1);
    nvec = n-mmax*tau; 
    if nvec-theiler < 2
    xM = NaN*ones(nvec,mmax);
    for i=1:mmax
        xM(:,i) = xV(1+(i-1)*tau:nvec+(i-1)*tau);
    % If theiler window define the components of the vector of all
    % distances that have to be removed.
    if theiler>0
        targV = [0:nvec-theiler-1]*nvec-isumV(1:end-mmax*tau);
        ntarg = length(targV);
        omitV = NaN*ones((nvec-theiler)*theiler,1);
        for ith=1:theiler
            omitV(([1:ntarg]-1)*theiler+ith) = targV+ith;
        omitV = [omitV' omitV(end)+1:(nvec*(nvec-1)/2)];
    for im=1:length(mV)
        % To handle large time series use single precision.
        if nvec>10000
            distV = single(pdist(xM(:,1:im)));
            distV = pdist(xM(:,1:im));
        % Remove components of distance matrix that regard temporally close
        % points
        if theiler>0
        % Count the distances that lie in each bin formed by the grid of
        % distances (given by 'resol')
        countV =histc(distV,rV);
        densityM(:,im) = countV/length(distV);
    % The log10 of distances corresponding to each bin
    log10rV = (resol-[1:resol]')*logdmin/(resol-1);
    cM = cumsum(densityM(1:end,:)); 
    cM(cM==0)=1; % To avoid log10(0)
    log10cM = log10(cM);
    log10rV = log10rV(2:end);
    log10cM = log10cM(2:end,:);
    rV = 10.^log10rV(1+smoothorder:resol-1-smoothorder);
    for im=1:nm
        m = mV(im);
        % Smooth derivative of the log10 of the correlation sum with 
        % respect to the log10 of radii.  
        % Compute the correlation dimension only if enough radii are
        % given to perform the smooth derivation
        poidimV = NaN*ones(resol-1-2*smoothorder,1);
        for i=1+smoothorder:resol-1-smoothorder
            A = [log10rV(i-smoothorder:i+smoothorder)';ones(1,2*smoothorder+1)]';
            b = log10cM(i-smoothorder:i+smoothorder,im);
            x = A \ b;
            poidimV(i-smoothorder) = x(1);
        for is=1:ns
            s = sV(is);
            % Compute the correlation dimension only if enough radii
            % are given to form scaling intervals according to 'e'.
            if rV(end) / rV(1) > s
                interval = length(find(rV<=rV(1)*s));
                nusd = m+1; % large starting values for mean and SD
                numean = NaN;
                for k=1:length(rV)-interval
                    nusdi= std(poidimV(k:k+interval));
                    numeani = mean(poidimV(k:k+interval));
                    if nusdi < nusd & numeani>0.7 
                        numean = numeani;
                        nusd = nusdi;
                nuT(itau,im,is) = numean;

Contact us