@ This is a Gauss code for the Computation of a Single PE point estimate of a univariate time series and associated Relative Frequencies @ @ This code accompanies the Aptech blog Permutation entropy which was published on 12/14/2018 @ @ See: https://www.aptech.com/blog/permutation-entropy/ @ @ This Gauss program requires a working copy of GAUSS 18+ @ @ Suggested Citation: Clower, E. and M. Henry, 2018. "petropy: Gauss module to compute a Single Permutation Entropy point estimate," @ NEW; CLS; /* ** Format: H = petropy(x, D, tau, method, accu) ** ** Input: ** x Column vector, univariate time series data ** ** D Scalar, embedding dimension ** ** tau Scalar, embedding time delay that determines the time separation between time series values ** ** 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 ** ** Outputs ** H Permutation Entropy ** Hnorm Normalized Permutation Entropy ** refreq Relative Frequencies with its plot */ // one-dimensional time series toy example x = {4, 7, 9, 10, 6, 11, 3}; x = packr(x); //embedding dimension D = 3; //embedding time delay tau = 1; // PE print; H = pentropy(x, D, tau, "noise"); print "Permutation Entropy point estimate =" H[1]; print; // RELATIVE FREQUENCIES rem = {1}; RelFreq = delRow(H, rem); print "Relative Frequencies" RelFreq; numrows = rows(RelFreq); xx = seqa(1, 1, numrows); plotXY(xx, RelFreq); #include dynargs.dec 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,Hnorm,factor; n_dynargs = COUNT_DYNARGS; //Set accuracy of the values in X if n_dynargs < 2; accu = 4; // default else; accu = sysstate(GET_ONE_DYNARG, 2); endif; //Set method if n_dynargs < 1; method = "order"; // default else; method = sysstate(GET_ONE_DYNARG, 1); endif; //Get length of x 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); // build n x (M-tau(n))-matrix from shifted values in x 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); // Column vector 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)'; 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; // INDEX ORDERS of smallest to largest sort_ind_mat = sortmatind(shift_mat); ind_mat = zeros(rows(sort_ind_mat), cols(sort_ind_mat)); // Filling the zero matrix with RANK ORDERS for ii(1, cols(ind_mat), 1); ind_mat[sort_ind_mat[., ii], ii] = seqa(1, 1, n); endfor; endif; // assign unique number to each pattern (base-n number system) 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); // Number of admissible Permutation Patterns permpat_num = diff(tmp); // RELATIVE FREQUENCIES of PERMUTATIONS p[i] refreq = permpat_num/sumc(permpat_num); // PE of order D H = -sumc(refreq .* log2(refreq)); factor = n!; Hnorm = (1/log2(factor)) * H; retp(Hnorm|refreq); 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;