c(1,1,1,2,2,2,3,3,3)

then

entropy.count

calculates the Shannon entropy using empirical/maximum likelihood probabilities for all unique symbols 1,2,3.

entropy.count <- function(entry) {
	counts <- lapply(split(entry, as.factor(entry)), length)
	counts <- unlist(counts)
	ps <- counts / sum(counts)
	entropy(ps)
}

The code is now also available through the "Papers / Code" tab above.

For python there is pyentropy

For R there is the entropy package by Hausser and Strimmer.

If you know of others, please let me know.

Also, I will soon add an additional page to the blog where you can additionally download all the source code presented in the code.

There's a vast number of papers and techniques on reducing the bias and variance on entropy estimates and I decided to write a few posts about this, with the aim to find the best entropy estimators for our (local) uncertainty measure. With a suitable entropy estimator we will be able to analyse local predictabilities conditioned on larger number of previous symbols with higher significance.

The estimator we used so far is called "plug-in" or maximum likelihood estimator and is defined as

where , so the number of occurrences of the word x in the whole space. It is well known that the MLE estimator is negatively biased. What does that mean?

Bias is the difference between the expected value of the estimator and the true value, i.e.

Thus it is a measure of how accurate the estimator is. The negative bias of MLE means that it underestimates the true entropy of a system. Let us check  if we can find this behaviour in data -- theory is nice but data is what counts.

We need to set up a system where we can easily check MLE with the true entropy, so it's easiest if we just sample from a uniform distribution with, e.g. 256 different observations, thus the true probability is for each observation; it easily follows that the entropy is then just .

The function biassim below calculates the bias through steps numbers of samples with replacement from a pool of 100'000 uniformly distributed integers between 0 and 256; it then calculates the entropy of the sample. Also notice here that the number of samples and number of possible words are equal, which is obviously  not the case in general. Finally, we subtract the true entropy value trueEV from the estimates.

biassim <- function(n, m, steps) {
	trueEV <- log2(m)
	uniformM<-round(runif(100000, min=0, max=m))
	f <- function(i) {
		sampleN <- sample(uniformM, n, replace = TRUE)
		probN <- sampleN / sum(sampleN)
                entropy(probN[probN > 0])
	}
	accu <- sapply(1:steps, f)
	return(accu - trueEV)
}

We apply the function above as follows:

> r <- biassim(256,256,2000)
> hist(r,breaks=100,xlim=0:-1,main="Bias and variance of MLE")

The histogram should like similar to the plot below.

It shows the negative bias very nicely, with the red vertical line being the true value to be estimated.

If our naive MLE estimator has this systematic underestimation then surely we can easily fix this issue by "shifting" the estimates to the positive side? Yes and no. The simplest bias correction is known under the Miller-Maddow correction and it suggests adding to the estimate. The only problem is that generally m, the number of possible words, is unknown, so it would have to be estimated as well.

In the next post we will apply the Miller-Maddow correction with a few other known estimators and analyse which one is most suitable for our purposes.

One important unanswered question so far is about the significance of the local uncertainties . Does a deviation from almost no order ( > 0.99) really mean something or is it due to imprecisions/undersampling of the empirical probabilities? As the original paper notices, the larger values we choose for n, i.e. the more previous trading days we consider to predict the next one, the more ngrams are possible and therefore the more samples we need to approximate the probabilities more or less accurately.

There's two ways to go:

• As in the original paper, use empirical probabilities and the basic plugin entropy estimator and restrict n to maximally 5, as their significance level K dictates (more to that below)
• Experiment with larger n including more sophisticated probability and enstropy estimators

We will do both. But for now I'll concentrate on the significance level K as introduced in the paper. A so called surrogate sequence of length n is generated out of the partitioned time series. These surrogates have the same mean and standard deviation as the original sequence, you could see it as a random shuffling of the sequence with some further rules. The local uncertainties from the surrogates are called . The significance level K is then calculates as:

where is the mean and is the standard deviation of the surrogate's answers. Intuitively, this is the difference between the actual uncertainties and mean randomised ones in units of the standard deviation of the surrogates. A K value larger than 2 would mean a confidence level greater than 95%. Due to the exponential nature however we'll just say, the bigger K the better. The reference to this test is from the Nature magazine 379 618 Characterisation of Low-Dimensional Dynamics in the Crayfish Caudal Photoreceptor.

It is not completely clear how this should get implemented. For now, I decided to create 5 surrogate sequences of the length n, and build up the definition above from that.

K <- function(seq, slideseq, ngram) {
	# generate 5 surrogate sequences out of seq
	sur1 <- surrogate(seq, length(ngram))
	sur2 <- surrogate(seq, length(ngram))
	sur3 <- surrogate(seq, length(ngram))
	sur4 <- surrogate(seq, length(ngram))
	p1 <- h_n_cond(sur1, ngram,3)
	p2 <- h_n_cond(sur2, ngram,3)
	p3 <- h_n_cond(sur3, ngram,3)
	p4 <- h_n_cond(sur4, ngram,3)

	s <- sd(c(p1,p2,p3,p4))
	#print(s)
	m <- mean(c(p1,p2,p3,p4))
	p <- h_n_cond(slideseq, ngram,3)
	return( abs( (p - m) / s) )
}

The length of the alphabet is fixed to 3 here in the code, but you can change that easily.

Now I'd like to calculate a sliding window of 4-grams to calculate local uncertainties of the last 200 DJI closes and their corresponding K values. The following function does just that

localallinfo <- function(seq, slideseq, lambda, limit=0:0) {
  len <- dim(slideseq)[1]
  range <- 1:len   if(length(limit) > 1) {
  	  range <- limit
  }
  # we store 2 columns for h_n_cond value and K, and dim(slideseq)[2] for length of word
  v <- array(dim=c(length(range), 2 + dim(slideseq)[2]))
  for(i in range) {
	v[i,] <- c(h_n_cond(slideseq, slideseq[i,], lambda),
                       K(seq, slideseq, slideseq[i,]), slideseq[i,])
  }
  return(v)
}

This will return a 2+n column array with the values for (local uncertainty, K value, ngram). Stored in variable x we have the return value of localallinfo. We can query all the results with high predictability:

> x[x[,1] < 0.932,]
          [,1]      [,2] [,3] [,4] [,5] [,6]
[1,] 0.9308985 30.705164    0    2    1    1
[2,] 0.9308985 13.576514    0    2    1    1
[3,] 0.9246361  6.738199    1    1    2    1
[4,] 0.9308985  7.542447    0    2    1    1
[5,] 0.9246361 13.721803    1    1    2    1
[6,] 0.9245277  9.663806    1    1    1    1

Very satisfying to see this develop so nicely. So what we have here are the patterns, their local uncertainty and the K values, all bigger than quite a bit bigger than 2. We can read this as follows for entry 1:

The ngram pattern 0 2 1 1 has a local uncertainty of 93% with a significance level K of 30.7

Thanks to the R package ggplot2, we generate the graph below with the commands:

> z<-data.frame(A=x[,1],B=1:200,K=x[,2])
> d<-qplot(B,A,data=z,colour=K,xlab="Time",ylab="Local Uncertainty")
> d + scale_colour_gradient(limits=c(1, 15), low="yellow",high="red")

This is the same data as the plot in the previous post, just now color coded with the significance level. We can see the trend that lower uncertainties have higher K values, just like described in the original paper.

In the next post I'll try to play a bit more with graphical representations, K levels, and then will move on to different entropy estimators and applying our code to more time series.

First we assume that already have our real returns data partitioned into symbols so is 3. Thus our time series is just a vector of values 0 1 2.

Next, all our functions will consider trajectories of that original vector. I will implement this as a sliding window of length n. So if our sequence is 012020120 the function slide will create the array 012, 120, 202, 020, 201, 012, 120 out of it.

slide <- function(seq,windowsize) {
	steps <- length(seq)-windowsize
	start <- 1
	stop <- windowsize
	accu <- array(0,dim=c(steps,windowsize))
	for(i in 1:(steps)) {
		#print(seq[start:stop])
		accu[i,] <- seq[start:stop]
		start <- start+1
		stop <- start+windowsize-1
	}
	return(accu)
}

The probability of a trajectory of length n is . This is simply modeled by the function prob_n by the empirical distribution of the number of times the ngram is found, divided by the total number of windows

prob_n <- function(seq, ngram) {
 l <- dim(seq)
 count <- 0
 for(j in 1:l[1]) {
	if(identical(seq[j,], ngram)) { count<-count+1 }
 }
 return(count/l[1])
}

The actual entropy is the average uncertainty in an ngram given the complete sequence. So if a specific ngram has probability then its weighted self-information is . This quantity is added up for all possible ngram sequences. This is implemented in the following function

H_n <- function(seq,n,lambda) {
	nseq <- slide(seq,n)
	# get the unique ngrams
	uniques <- unique(nseq)
	dim_u <- dim(uniques)
	accu <- 0
	# for every unique ngram, get its prob_n and add to its entropy
	for(i in 1:dim_u[1]) {
	  p <- prob_n(nseq, uniques[i,])
	  #accu <- accu + p
	  accu <- accu + (-p *log(p)/log(lambda))
	}
	return(accu)
}

The average uncertainty of predicting the next symbol , given that you've already seen n previous ones is just the difference of as follows

h_n <- function(seq,n,lambda) {
	return(H_n(seq,n+1,lambda) - H_n(seq,n,lambda))
}

The more interesting bit is how to calculate the conditional probabilities in , . Say the ngram we are looking at is the sequence 0120, is the uncertainty in predicting the next symbol after observing 0120. This is implemented in a way where we only consider the sequences of length 5 with an unknown fifth symbol 0120x. The code h_n_cond finds the empirical probability distribution for x in the restricted sample space of 0120x. Once the probability distribution is accumulated in the variable probdist, we simply calculate the based entropy of it to arrive at .

h_n_cond <- function(seq, ngram, lambda) {
  nextcounts <- array(data=0, dim=c(lambda,1))
  l <- dim(seq)
  skip <- length(ngram)
  #count <- 0
  for(j in 1:l[1]) {
  	if(identical(seq[j,], ngram)) {
  	   #count <- count + 1
  	   # keep count for each possible next step element
  	   if(j < l[1]-(skip-1)) {
  		nextcounts[seq[j+skip,1]+1] <- nextcounts[seq[j+skip,1]+1] + 1
  	   }
  	}
  }
  probdist <- (nextcounts/sum(nextcounts))
  entropy(probdist, lambda)
}

Testing

Next up, we need some real data to test our functions on! I downloaded all of DJI historical data which was available from Yahoo. The data is stored in the variable dji and I will only consider dji$Close prices for now. I generate the daily logarithmic returns which will be partitioned into the discrete symbols 0, 1 and 2. We will use the same partitioning as in the paper. In R this is quite straightforward as:

> djipart <- ifelse(djiret < -0.0025, 0, ifelse(djiret > 0.0034, 2, 1))

The authors justified the asymmetric partitioning due to the slightly positive mean logarithmic prices. I rather see a negative mean of -0.0001835142 in my data but I want to reproduce the numbers of the paper as well as possible so I use theirs.

Before running any functions on the real data I want to test it on uniformly distributed random data

> sequnif <-round(runif(10000, min=0,max=100)) %% 3

Obviously, I expect any uncertainties calculated as equally likely .. there shouldn't be any higher predictability by considering longer trajectories

> h_n(sequnif, 3, 3)
[1] 0.997702
> h_n(sequnif, 4, 3)
[1] 0.9928539
> sequnif4 <- slide(sequnif,4)
> h_n_cond(sequnif4, c(1,0,1,0), 3)
[1] 0.9980187
> h_n_cond(sequnif4, c(1,0,1,1), 3)
[1] 0.9982527
> h_n_cond(sequnif4, c(0,1,1,1), 3)
[1] 0.9974415

Our functions seem to produce the expected almost total uncertainty answers. The conditional entropies for n from 1-6 is plotted in the picture below. Here are some numerical values:

> h_n(djipart, 3, 3)
[1] 0.9841248
> h_n(djipart, 4, 3)
[1] 0.9788154
> h_n(djipart, 5, 3)
[1] 0.9679176
> h_n(djipart, 6, 3)
[1] 0.9429579

We observe a nice gradual decrease of uncertainty for longer trajectories.

This graph pretty much reflects the shape of Fig. 2 in the original paper. Again, the graph reflects the average uncertainty of predicting the next symbol for ngrams of different lengths (with values from the x axis).

The local uncertainty for specific ngrams can be queried with h_n_cond as follows:

> h_n_cond(dji4, c(0,2,1,2), 3)
[1] 0.9924063
> h_n_cond(dji4, c(1,1,1,1), 3)
[1] 0.9245277

The next state after observing pattern 1111 is more predictable than the pattern 0212 by 6.7%!

The next graph I will generate is the local uncertainty of length 4 of the last 200 trading day closes using a rolling window as follows

len <- dim(dji4)[1]
v4 <- vector()
for(i in 1:200) {
   h <- h_n_cond(dji4, dji4[i,], 3)
   v4 <- c(v4,h)
}

> plot(rev(v4), type='l', axes=F, xlab="DJI close", ylab="local uncertainty")
> axis(1,at=pretty(c(1,200)),labels=rev(dji$Close[pretty(c(0,201))]))
> axis(2)

Which resembles the graph from Fig. 1 from the original paper. Given the original definitions we should have that the mean of that local uncertainty sequence v4 should be close to the average uncertainty of length 4:

> mean(v4)
[1] 0.9764447
> h_n(djipart, 4, 3)
[1] 0.9788154

The minimal value in this graph is at 0.9245.

What we haven't considered so far is how significant the predictions are. The MLE estimator of entropy we are using here is biased and the sample error grows rapidly with increasing ngram length. The paper is using something they call surrogate sequences to calculate its significance. I don't really know about that approach but I will look into it and write about it in a future article in this series of local uncertainties. One approach I could imagine is using different entropy estimators with bias correction.

Hope you enjoyed this article so far, please leave a comment if you feel like it.

The paper presents an excellent application of information theory to time series analysis. The idea is simple: is it possible to find sub-trajectories in financial time series (here the daily returns of some indices or stock) where a "local order" exists with higher than average predictability.

I won't explain the paper in full, so please have a look at the pdf above for notation and details. However I will describe the most important concepts below. We consider one-dimensional, discretely partitioned time series. The authors use Shannon entropy H as basic tool to measure uncertainty or predictability of the probability distribution described by the time series. For a certain trajectory of length n the uncertainty of predicting the next state is the difference in Shannon entropies for trajectories of length n+1 and n:

Let be the total length of the time series, let be the length of the alphabet and is a specific subtrajectory of length .

Further, is the probability to find a subtrajectory with the letters in the complete time series. Then the entropy is simply

Notice here that we take the logarithm of base , so the quantities are not measured in the more conventional bits. As described above from the conditional entropy we define the average predictability as

Clearly, since we consider units, the maximum uncertainty is . Also we can deduce that which is intuitively understandable as: the more I know the less uncertainty there can be, unless of course every observation was completely independent, where the equality would apply.

Next, the paper describes the "local uncertainty" measure as the uncertainty of the next symbol given the observation of specific symbols . The authors write it as

Notice that this almost coincides with the conditional entropy except that it leaves out the weighting factor of , thus it's not an averaging measure but a local one. The function will be the main measure we apply to our time series.

Partitioning

As a next step we need to partition the real valued data of the time series into discrete symbols. First we take the daily logarithmic prices changes of some stock :

The returns now need to be partitioned into symbols of an alphabet with length . The optimal choice of partitioning is a whole other question, but for now we just assume that it doesn't matter. The paper uses for and partitions Dow Jones daily returns with the following thresholds:

This is obviously different for every stock we look at and needs to be chosen depending on its basic statistics (skew, mean).

A randomly chosen snippet of a trajectory then might look like the string 0112012220201.

This is all we need to know for our implementation, which will be described in the next post of this series.

So here's the fairly obvious way on how to calculate Shannon's entropy in R using Reduce:

> fentropy <- function(x,y) { x + (-y * log2(y)) }
> Reduce(fentropy, c(0.5,0.5), 0)
[1] 1
> Reduce(fentropy, c(0.25,0.25,0.25,0.25), 0)
[1] 2

First for the binary case with answer 1, and then for four values uniformly distributed.

Last but not least, we could also write an entropy function the "R way" which uses its nice functions which work over vectors:

entropy <- function(ps) {
     H = -sum(ifelse(ps>0, ps * log2(ps), 0))
     return(H)
}

The current list is:

Untertainty analysis in financial markets: can entropy be a solution?, Andreia Dionísio, Rui Menezes and Diana A. Mendes (pdf)

Forecasting Foreign Exchange Market Movements via Entropy Coding, Arman Glodjo, Campbell R. Harvey (pdf)

Local order, entropy and predictability of financial time series, L. Molgedey and W. Ebeling (pdf)

These three papers all use Shannon Entropy in place of more traditional statistical measures. What is interesting however is that all of them apply entropy in different ways.
The first paper by Dionisio compares entropy as measure of uncertainty with variance/standard deviation in portfolio management.
The second paper by Glodjo applies techniques from coding theory (the original and most successful application of IT) to forecasting high frequency time series. Also it provides good arguments for using IT in finance.
The last paper by Molgedey is using conditional entropy directly on returns time series to quantify "local order" in highly stochastic time series. A local order would be a point in time where the next step is more predictable than average.

I will have a look at some of those techniques in more detail and might implement some of it to see if I can replicate the authors results.

A left fold [foldleft f accu l] applies the head of the list l to the function f together with the accumulator variable accu. The result is the new accumulator which is used in the next recursive call together with the tail of the list.

A left (or right) fold is easily implemented in R as follows:

foldleft <- function(f,accu,l) {
	if(length(l)==0) {
		accu
	} else {
		head <- l[1];
		tail <- l[-1];
		foldleft(f, (f(accu, head)) , tail)
	}
}

To see how it works, we could apply it to calculate the variance of a fair die. Remember the variance is just where is the mean, which is implemented in the following function f:

mean<-sum(1:6)/6

f <- function(accu,i) {
	accu+(1/6 * (i-mean)^2)
}

foldleft(f,0,c(1,2,3,4,5,6))

where our last call to foldleft evaluates to 2.91667.