pentropy_rw_f.gss

new;
cls;
et  =  hsec  ; 

// DJIA Data
// 1901 - 2016
row_range={1,31519};
  
x = csvReadM("DJIA_1901_to_2016_red.csv",row_range,2);
x = packr(x);

y = ln(x[2:rows(x)]) - ln(x[1:rows(x)-1]);
x=y;

numr = rows(x);
print "Number of Obs.: " numr;
print;
 
numprom = 749;
// number of segments with a day shift in the matrix x_seg below
numb = numr - numprom + 1;  
print "Number of segments: " numb;

// Rolling Window
year=3;      
win = 750;

// Segmenting the data
x_seg = segmentEntropyData(x, win);

proc(1) = segmentEntropyData(data, win);
   local mat;
    mat = zeros(win,rows(data)-win+1);
    for i(1, rows(data)-win+1, 1);
        mat[.,i] = data[i:i+win-1];
    endfor;
        
retp(mat);
endp;

//embedding dimension
n = 5;
print "Permutation order:" n;
print;

//Time delay
tau = 1;
print "Time delay:" tau;
print;
       
H_mat = tvPentropy(x_seg, n, tau, "noise");     
H_matcolv=H_mat';

print "Min H: " minc(H_matcolv);
print;
print "Max H: " maxc(H_matcolv);
print;
print "Mean H: " meanc(H_matcolv);
e = {.025,.5,.975};
y = quantile(H_matcolv,e);
print "medians" ;
print y[2,.];
print;
print "95 percentiles";
print y[1,.];
print y[3,.];
print;

// Permutation entropy row vector
print "Permutation Entropy Statistic =" H_mat[1,.];     
print;

//Plot control structure
struct plotControl myplot;
myplot = plotGetDefaults("xy");

plotSetTitle(&myPlot, "Time-varying Permutation Entropy", "Arial", 18);
plotSetXLabel(&myPlot, "Obs.", "Arial", 12);
plotSetYLabel(&myPlot, "Normalized Permutation Entropy", "Arial", 12);
plotSetGrid(&myPlot, "off");

// Time varying entropy statistic
plotXY(myPlot, seqa(1,1, cols(H_mat)), H_mat[1,.]');

//Create time vector
dtstart = 20070102000000;
dtend = 20091231000000;

//Get dates
dates = getDates(dtstart, dtend);

proc(1) = tvPentropy(x_mat, n, tau, ...);
    local h_mat, tmp;
    
    h_mat = {};
    
    for i(1, cols(x_mat), 1);
    
        tmp = pentropy(x_mat[.,i], n, tau, ...);
        h_mat = h_mat~tmp;
    endfor;
    
    retp(h_mat);
endp;


#include dynargs.dec
/*
**
** Purpose: The petropy function calculates the permutation entropy of data series.
**
** Format:  H = PETROPY(X, N, TAU, METHOD, ACCU)
**
** Input:
**  X       Vector, data
**
**  N       Scalar, permutation order
**
**  tau     Scalar, time lags
**
**  method  String treatment of equal values
**          'noise' - add small noise
**          'equal' - allow same rank for equal values
**          'order' - consider order of appearance (first occurence --> lower
**            rank)
**
**  accu    Scalar, parameter describing accuracy of the values in X
**          by number of decimal places
**
**  Output
**  H      Permutation entropy
**
**
** References:
**
** Bandt, C.; Pompe, B. Permutation Entropy: A Natural Complexity
** Measure for Time Series. Phys. Rev. Lett. 88 (2002) 17, 174102
**
** Riedl, M.; M¸ller, A.; Wessel, N.: Practical considerations of
** permutation entropy. The European Physical Journal Special Topics
** 222 (2013) 2, 249ñ262
**
*/
proc (1) = pentropy(x, n, tau, ...);
    local n_dynargs, accu, method, M, equal, tmp, shift_map,
    k, j, shift_mat_ind, tmp2, tmp3, shift_mat, ind_mat, ind_vec, sort_ind_mat, ia, ic,
    tmp_sort, permpat_num, refreq, size,tmp0, denom, H, Hnorm,factor;
    
    n_dynargs = COUNT_DYNARGS;
    
    if n_dynargs < 2;
        accu = 4;           
    else;
        accu = sysstate(GET_ONE_DYNARG, 2);
    endif;
    
    if n_dynargs < 1;
        method = "order";   
    else;
        method = sysstate(GET_ONE_DYNARG, 1);
    endif;
    
    M = rows(x);
    
    equal = 0;
    
    if n*log(n)>15;
            errorlogat "Permutation dimension too high";
        end;
    endif;
    
    if (rows(tau)) > 1 and (rows(tau) != n-1);
        errorlogat "Time lag vector has to have n-1 entries";
        end;
    endif;
        if ((n-1)*minc(tau) >=M) or maxc(tau) >= M;
        errorlogat  "Too few data points for desired dimension and lags";
    endif;
    
    if lower(method) == "noise";
        print "Method: Add small noise.";
        x = x + rndn(M, 1)*10^(-accu-1);
    elseif lower(method) == "equal";
        print "Method: Allow equal ranks.";
        equal = 1;
    elseif lower(method) == "order";
        print "Method: Consider order of occurrence.";
    else;
        errorlogat "Unknown method. Default method 'order' used";
    endif;
    print;
    
    if rows(tau) > 1;
        tau = reshape(tau, rows(tau), 1);
        tmp = 0|tau;
        tau = sortc(tmp);
        
        shift_mat = zeros(n,  M-tau[n]);
        
        for ii(1, n, 1);
            k = tau[ii] + 1;
            j = M - tau[n] + tau[ii];
            shift_mat[ii, .] = x[k:j];
        endfor;
        
    else;
        
        shift_mat_ind = zeros((M-(n-1)*tau)*n, 1);
        tmp0 = seqa(1, tau, n);   
        shift_mat_ind[1:n] = tmp0;
        for i(1, (M-(n-1)*tau)-1, 1);   
            shift_mat_ind[(i*n+1):((i+1)*n)] = tmp0 + i;
        endfor;
        
        shift_mat = reshape(x[shift_mat_ind], (M-(n-1)*tau), n)';   // D x (M-D+1) matrix

    endif;
    
    if equal;
        ind_mat = zeros(size(shift_mat));
        for ii(1, cols(ind_mat), 1);
            tmp = unique(shift_mat[., ii]);
            ind_mat[., ii] = indnv(shift_mat[., ii], tmp);
           
        endfor;
    else;
        
        sort_ind_mat = sortmatind(shift_mat);                       
        
        ind_mat = zeros(rows(sort_ind_mat), cols(sort_ind_mat));    
        
        for ii(1, cols(ind_mat), 1);    
            ind_mat[sort_ind_mat[., ii], ii] = seqa(1, 1, n);       
        endfor;
    endif; 
    
    ind_vec = n.^(seqa(0, 1, n))' * (ind_mat - 1);    
    
    tmp_sort = sortc(ind_vec', 1);   
    
    tmp2 = unique(tmp_sort);   

    ia = indnv(tmp2, tmp_sort);     
    
    tmp = ia|(cols(ind_vec) + 1);     
    
    permpat_num = diff(tmp);                 
    refreq = permpat_num/sumc(permpat_num); 
    H = -sumc(refreq .* log2(refreq));
    factor = n!;
    print factor;
    Hnorm = (1/log2(factor)) * H;
    retp(Hnorm);
endp;

proc (1) = log2(x);
    retp(log(x)/log(2));
endp;

proc (1) = diff(x);
    local tmp;
    tmp = x[2:rows(x)] - x[1:rows(x)-1];
    retp(tmp);
endp;

proc (1) = sortmatind(x_mat);
    local sort_ind_mat;
    
    sort_ind_mat = zeros(rows(x_mat), cols(x_mat));
    
    for i(1, cols(x_mat), 1);
        sort_ind_mat[.,i] = sortind(x_mat[.,i]);
    endfor;
    retp(sort_ind_mat);
endp;

proc (1) = delRow(x, remove);
   local mask, xout;
   mask = zeros(rows(x), 1);
   mask[remove] = ones(rows(remove), 1);
   xout = delif(x, mask);
   retp(xout);
endp;

proc (1) = getDateRange(X, range);
    local ndigits, idx, X_range;
    ndigits = floor(log(range)) + 1;
    range = range .* 10^(14-ndigits);
    idx = sumc((x[.,1] .< range[1]) ~ (x[.,1] .< range[2])) + 1;
    X_range = X[idx[1]:idx[2],.];
    
    retp(X_range);
endp;

proc(1) = getDates(dtstart,dtend);
    local date_vec, tmp1;   
    
    date_vec = {};
    do while dtstart<=dtend;
        tmp1 = dtstart;
        dtstart = getNextTradingDay(dtstart);
        date_vec = date_vec|tmp1;
        
    endo;
    
    retp(date_vec);
endp;

print; 

print;
(hsec-et)./100 "seconds";