% Cascade model

%Reference:
% Flores, Claudia, Multiplicative cascade models for rain in
% hydro-meteorological disasters risk management

% Begin
m=4; % This is the mass to be distributed along the support.
n=8; % This is the number of levels of the cascade that we want to develop.
T=1; % This is the length of the support. The cascade has support in the interval [0,T].
branchn=2^n; % This is the branch number. In the last level (n) we will have 2n branches.
% Now we generate a matrix with complementary weights (w).
w=zeros(n,branchn);
for i=1:1:n; % rows
    for j=1:2:2^i %columns;
        w(i,j)=rand(1,1); % The weights are uniformly distributed in [0,1].
        w(i,j+1)=1-w(i,j); % We have complementary weights.
    end;
end;
% In the matrix M, we store the weights in a convenient order.
M = ones(n+1,branchn);
% The first row of M will have ones, because the mass (m) is uniformly distributed along the
% support at the first level of the cascade,
% The cells of the other rows will store elements of the matrix of weights (w).
for i = 2:1:n+1 % i-1 is the corresponding row in w
    repeat = 2^(n-i+1);
    for k = 0:1:branchn/repeat-1;
        for j = k*repeat+1:1:(k+1)*repeat;
            if j <= 2^k*repeat;
                temp = k+1;
            else temp = k+2;
            end;
            M(i,j) = w(i-1,temp);
        end;
    end;
end;
% We calculate the values for every branch:
cascade = cumprod(M)*m;
% Every row of the matrix 'cascade' corresponds to a level of the cascade.
% For the graphics:
N = zeros(n+1,branchn);
for i=1:1:n+1; % We should consider the different scales for the graphic.
    N(i,:)=cascade(i,:)*2^(i-1)/T; % . i-1 is the level of the cascade.
    % The calculated mass for every sub-interval of length T/2i-1 is distributed along it.
end;
graf=N';
t=T/branchn:T/branchn:T; % The matrix stores the support divided in the finest scale.
% This is to obtain a graphic of the simulated binary multiplicative cascade
figure(1);
for i=1:1:n+1;
    subplot(3,3,i); % n+1=9 graphics are generated and presented in a 3x3 array.
    plot(t,graf(:,i));
    axis([1/branchn,T,0,max(graf(:))]);
end;