% logistic map
a = [3.5:.001:4];
x(:,1:length(a)) = .6;
for i = 1:2020;
    for j = 1:length(a)
        x(i+1,j) = a(j) * x(i,j) * (1 - x(i,j));
    end
end
x(1:21,:) = [];




% make symbol sequence (consider only one time series)
i = 450; % for a(i)
s = x(:,i) > .7;


% probalities of symbols
p_1 = sum(s==1)/length(s);
p_0 = sum(s==0)/length(s);

% joint probabilities
for t = 0:40
    p_11(t+1) = sum(s(1:end-t) == 1 & s(1+t:end) == 1)/length(s);
    p_10(t+1) = sum(s(1:end-t) == 1 & s(1+t:end) == 0)/length(s);
    p_01(t+1) = sum(s(1:end-t) == 0 & s(1+t:end) == 1)/length(s);
    p_00(t+1) = sum(s(1:end-t) == 0 & s(1+t:end) == 0)/length(s);
end

% conditional probabilities
p_1_1 = p_11./p_1;
p_1_0 = p_10./p_0;
p_0_1 = p_01./p_1;
p_0_0 = p_00./p_0;



% mutual information
I(i,:) = .5 * (p_1_1 .* log(p_1_1/(p_1*p_1)) + ...
p_1_0 .* log(p_1_0/(p_1*p_0)) + ...
p_0_1 .* log(p_0_1/(p_0*p_1)) + ...
p_0_0 .* log(p_0_0/(p_0*p_0)));





% word statistics
symbols = 2; % number of symbols
l = 3; % word length

N = symbols^l; % number of all possible words
w = dec2base(0:N-1,symbols); % set of all words of length l


ss = num2str(s)'; % symbol seqz
uence as a string (because w is a string matrix)

% histogram of word occurrences
for j = 1:length(w)
    wn(j) = numel(strfind(ss, w(j,:)));
end

% probability of word occurrences
pn = wn/nansum(wn(:));





% example Rössler (ode integration, external file roessler.m needed, see below)
% periodic oscillation
a=.2; b=.7; c=3.6; 
[t1,x1] = ode45('roessler', [0:.1:600], [16 8 4], [], [a b c]);
% chaotic behaviour
a=.2; b=.2; c=5.7; 
[t2,x2] = ode45('roessler', [0:.1:600], [16 8 4], [], [a b c]);






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% save the following line as an external file with the name roessler.m

function dy=roessler(t, y, dummy, p)
% periodic
% a=.2; b=.7; c=3.6;
% chaos (standard)
% a=.2; b=.2; c=5.7;
% a=.1; b=.1; c=9;
% funnel
% a=.2925; b=.1; c=8.5;

dy = zeros(3,1);

if nargin < 3 | isempty(p)
   a = .2;
   b = .2;
   c = 5.7;
else
   a = p(1);
   b = p(2);
   c = p(3);
end

dy(1) = -y(2) - y(3);
dy(2) = y(1) + a * y(2);
dy(3) = b + y(3) * (y(1) - c);