Quantitative thoughts

Few weeks ago GitHub announced, that its timeline data is available on bigquery for analysis. Moreover, it offers prizes for the best visualization of the data. Despite my art skills and minimal chances to win beauty contest, I decided to crunch GitHub data and run data analysis.

After initial trial of bigquery service, I found hard to know, what price, if any, I’m going to pay for the service. Hence, I pulled the data (6.5 GB) from bigquery on my machine and further I used my machine for analysis. Bash scripts have been used to clean up and extract necessary data, R for data analysis and visualization and C++ for text extraction.

GitHub dataset is one table, where each row consist of information about repository (i.e. path, date of creation, name, description, programming language, number of forks/watchers and etc.) and action, which was done by user (i.e. username, location, timestamp and etc.).

As a result, we can check how GitHub users actions are spread over time during the day. The X axis on the graph below is labeled with the hours of the day (GMT) and the Y axis represent median values of the actions for each hour. From it, we can make a deduction, that highest load for GitHub can be expected between 15:00 and 17:00 GMT and lowest to be expected between 05:00 and 07:00 GMT. The color of the line indicates how busy was the day based on quantiles: green are calm days (20% of days), blue – normal days (50% quantile) and red are busy days (80% quantile). I should to mention, that auto-correlation or serial correlation is high (70% for following hour), which means, that busy hours tend to be followed by busy hours and calm hours tend to be followed by calm hours. Moreover, busy days tend happen after busy days.

Second graph below shows median of actions divided by weekdays. There is not big surprise – weekends are more slow than weekdays, nevertheless the programmers are slightly less productive on Mondays and Fridays.

The analysis of creation of new repository shows, that the pattern of busy or calm hours remains over the years. This can be attributed to the fact, that majority of the users comes from North America and Europe.
Another hypothesis can be drawn from this information, that number of creation of the new repositories grow exponentially. However, I mind you, that the graph below is biased – most likely, GitHub users update recent projects, consequently more recent projects appeared on timeline. Even though, 2009-2011 years show exponential grow.
The X axis of the graph below is labeled with the hour of the day, the Y axis – log of median values of new repositories.

Following graph shows the number of forks per project (the X axis, log scale) versus number of watchers (the Y axis, log scale). As expected, there is linear correlation between forks and watchers. Even so there is something interesting about outliers, which are below bottom line – the projects, where number of watchers is low, but number of forks is high. These are anomalies and worth to check.

The next thing to do is to look at the repository description. Let’s group the repositories by programming language and count most dominant words in the description. The graph below has C++ word cloud on the left and Java – right . C++ projects are about library, game, simple(?), engine, Arduino. Java is dominated by android, plugin, server, minecraft, spring, maven.

Ruby (left) vs Python(right ):

“Surprise”, “surprise” – R projects (left) are largely about data analysis, however “machine” word, which corresponds to Machine learning is very tiny. Shell (right) is dominated by configuration, managing, git(?).

GitHub dataset includes location field. Unfortunately, the users can enter whatever they want – country, city or leave it empty. Nevertheless, I was able to extract good chunk of actions, where location field has meaningful value.  The video below shows country based users activity, where dark red corresponds to high activity and light red – minor. Only 30 most active countries are included, the rest are grey.
The same pattern persist over the days – activity in Asia increases around midnight, Europe wakes up around 8:00 or 9:00, where America starts around 15:00. Who said, that hackers and programmers work at night?

What else can be done with GitHub dataset? Most repositories have description field, which can be used to find similar projects by implementing tf-idf method. I tried that method and the results are satisfying.

Most of the graphs shown above are reproducible (except word clouds) and the code can be found on GitHub.

There are plenty of data on internet, however it is raw data. Think for a second about public surveillance cameras - useful to check the traffic on the route or busy place, but anything else? What if you want to know how many cars are on the route? How many car were yesterday at the same time? Given so many cars on the route, how much polluted air in the area?
While working on the road map for data dive event, I started to wonder, how feasible is to use data of public surveillance cameras. So I quickly built a pilot project and now I would like to share my experience.

First step – data acquisition. At beginning I was thinking to plug my smartphone somewhere and collect data of the busy route.  Nevertheless, I quickly found surveillance cameras in Vilnius and started to collect images. Run a search and I’m sure, that you will find them in your city:

Here is bash script, which I use to collect images:

#you need full path for crontab
cd /home/git/carCount/img
a=`date +%s`
b=${a}_4.jpg
wget -O $b -q "http://www.sviesoforai.lt/map/camera.aspx?size=full&image=K7742-1.jpg&rnd=0.15417794161476195"

Data preparation. After while you will have enough data to train your machine (for beginning more than 30 images should be O.K.).
How do we train the algorithm? The goal is to identify the cars in a given image. That means, that we have to provide the examples of positive images (clear image of the cars) and negative images (no car, parts of the car and etc.). Important note – we don’t feed whole image, but we cut a chosen image with sliding window (100×100 in my case). 4 examples of positive images:

Meanwhile, it is worth converting each image to portable grey format PGM. For this specific task, we can sacrifice information about the color of the car – it won’t improve prediction. Besides, PGM images can be loaded into R and easily transformed into matrix. Here is bash script, which converts jpg to pgm and slices each image:

#remove image duplicates
find . -maxdepth 1 -type f -exec md5sum {} \;  >test.txt
awk 'a[$1]++ {gsub(/^\*/,"",$2); print "rm ", $2}' test.txt |sh
rm test.txt
 
#convert jpg
 
if [ -d "out" ]; then
        rm -r out
fi
mkdir out
for k in $(ls *.jpg); do convert $k out/$k.pgm; done
 
cd out
mkdir slide
for filename in $(ls *.pgm);
 do 
 
w=`convert $filename -print "%w" /dev/null`
h=`convert $filename -print "%h" /dev/null`
let "ww= $w/100"
let "hh= $h/100"
for((y=150;y<=250;y+=50))
do
for((i=100;i<=400;i+=50))
do
echo "slide/$i.$filename"
let "h_slide=$i"
convert $filename -crop 100x100+$i+$y slide/$y.$i.$filename
done
done
done

Training, predicting, cross validation. Now is time to open R, load 100×100 images from “train/out/slide” directory and train the algorithm. Important note – each image is a matrix, however you have to feed a matrix of all images to learning algorithm (support vector machine in my case). What you have to do is to “unroll” each image matrix into a vector, get 1X10000 vector and build a new matrix, where each row is an image.
Once training is done, load unseen data from “crossval/out/slide” directory and check “result/” directory, where you will find  images of the cars. R script, which does all above:

setwd('/home/git/carCount/')
 
######read positives############
files=list.files('test/pos/')
pos=matrix(nrow=NROW(files),ncol=100*100)
 
for(i in 1:NROW(files))
{
  gray_file=read.pnm(paste('test/pos/',files[i],sep=''))
  pos[i,]=c(gray_file@grey)
}
outcome=vector(length=NROW(files))
outcome[which(outcome!=1)]=1
 
########read negatives#############
files=list.files('test/neg/')
neg=matrix(nrow=NROW(files),ncol=100*100)
 
for(i in 1:NROW(files))
{
  gray_file=read.pnm(paste('test/neg/',files[i],sep=''))
  neg[i,]=c(gray_file@grey)
}
tmp=vector(length=NROW(files))
tmp[which(tmp!=0)]=0
outcome=c(outcome,tmp)
forecast=svm(rbind(pos,neg),outcome)
cross_val=pos[84:90,]
pred=predict(forecast,cross_val,decision.values=TRUE)
 
##########################unseen data######################
files=list.files('crossval/out/slide/')
cross=matrix(nrow=NROW(files),ncol=100*100)
 
for(i in 1:NROW(files))
{
  gray_file=read.pnm(paste('crossval/out/slide/',files[i],sep=''))
  cross[i,]=c(gray_file@grey)
}
pred=predict(forest,cross,decision.values=TRUE)
 
###############copy positives into result directory###############
dir.create('result')
file.copy(paste('crossval/out/slide/',files[which(as.double(pred)>0.6)],sep=''),'result/')

Classified as positive by algorithm:

Classified as negative by algorithm:

Conclusion. It is truly amazing how well algorithm is able to separate wheat from the chaff without additional tuning. Mind you, my impression is biased after so many fails with financial data, which is noisy and good predictions are scarce.
Nevertheless, this project is far away for ideal – it doesn’t take into account weather condition, traffic jams, perspective view, movements of the camera and etc. But I leave this fun for data-dive event.

Fork the code: https://github.com/kafka399/carCount/

The first thing, what you have to do is to estimate how many users will be interested in local R group. I would say, that out of one million inhabitants you can expect 10-20 users. Based on this raw number, you can know, what challenges are waiting for you. If you expect 100 or more users, you have to think about the appropriate place to hold first meetup, how to manage so many people, what topics to present first time.
However, if you expect small community (as I did), your challenge is spreading the news about local R group. Get to know a few local users and ask them, what do they think about a meetup. The right place to find such users would be local university. Most likely the local university will be happy to provide a place for the first official meeting.

Fortunately for me, I met a powerful R user – @mpiktas, which is a lecturer at Vilnius University. With his help I was able to identify more R-infected users and get premises at university for presentation.

Next step, when you are sure, that you are not alone, is to choose a name for the group, build simple web site from scratch or by using a service like meetup. This won’t cost you a fortune, however I need to say, that you can apply for sponsorship. Revolution Analytics not only provides sponsorship, but as well manages a dictionary of R user group.

Once you have created virtual community you have to thinking about a meetup. As a founder, prepare an introduction about local R group, its future plans, your own usage of R. Here mine, which I used for VilniusR introduction:

What’s next? During first meetup you can outline future meetup, however I found fascinating, that many participants expressed an interest in datadive event. So, now we are in the process of organizing such event!

p.s. if you happen to be from Luxembourg and you are interested in local user group – let me know.

At work, the client requested, if existing search engine could accept singular and plural forms equally, e. g. “partner” and “partners” would lead to the same result.

The first option – stemming. In that case, search engine would use root of a word, e. g. “partn”. However, stemming has many weaknesses: two different words might have same root, a user can misspell the root of the word, except English and few others languages it is not that trivial to implement stemming.

Levenshtein distance comes as the second option. The algorithm is simple – you have two words and you calculate the difference between them. You can insert, delete or replace any character, but it will cost you. Let’s imagine, an user enters “Levenstin distances” into search engine and expects to find revalent information. However, he just made 2 errors by misspeling the author’s name and he used plural form of “distance”. If search engine accepts 3 errors – the user will get relevant information.

The challenge comes, when you have a dictionary of terms (e. g. more that 1 mil.) and you want to get similar terms based on Levenshtein distance. You can visit every entry in the dictionary (very costly) or you can push dictionary into the trie. Do you need a proof for the cost? There we go:

Red color indicates the performance of the search, when all terms are in the trie, green – simple dictionary.

Now we come to the second part of the post – why to bother and plot such graphs, if we could check few entries to determine average time and the winner? The reason is simple – we trust in God, all others must bring data. To say it differently – while profiling the code, you should be interested in average time AND variation. As you can see in the graph above, variation of the blue color is very small – it takes approximately the same time to scan whole dictionary. However, red has higher variation – the result can take for while or it can finish just at the beginning, but overall it works faster.
Now, imagine, that a programmer wants to define, which implementation A or B for volatile cache is much faster. Let’s assume, that big O notion is not going to help and she conducts 2 test for A and 2 for B. While running test A, cache size expands, while B – shrinks. As the result, B wins over A and she makes wrong choice. However, her colleague claims, that despite A has greater volatility, it is much faster and she tried with 500 queries! Whom should I trust?

I use this piece for code profiling:

1
2
3
4
5
6
7
8
9
10
11
12
13

 simple=read.table('simple.txt')
node=read.table('node.txt')
 
simple=cbind(simple,as.character(c('simple')))
colnames(simple)=c('time','type')
node=cbind(node,c('node'))
colnames(node)=c('time','type')
 
rez=data.frame(rbind(simple, node))
 
require(ggplot2)
 
ggplot(rez,aes(time,fill=type))+geom_density(alpha=0.6,size=1.3)+scale_x_log10()

The data, C++ code for Levenshtein distance and trie can be find on GitHub.

I found this source very useful: http://stevehanov.ca/blog/index.php?id=114

In the previous post I shared an example how to get high frequency data from IB broker (well, it is retail version of HFD – it has only best bid/ask and the trades). Now, once you saved some data – what should you do next?

Next logical step would be data sanity check and visualization. For example, while preparing R script for this post, I found, that IB data contains numerous duplicates in the quotes. Every time, when the trade happens, IB trading platform sends the price and the size of the trade bundled together. Additionally, it sends the size of the trade as separate quote as well and this completely mess up the data. So, data sanity check and visualization gave me a hint, that something is wrong with the data.

Today I want to show an example in R, which loads data from mongodb and plots some parts of the data. This should give you better intuition on collected data.

The plot shows bid (light blue) prices , ask(green) prices and the trades (red). The size of the red dot indicates volume of the trade.

The source code is shared on github and below:

#Author Dzidorius Martinaitis
#Date 2012-03-01
#Description 
 
require(rmongodb)
require(xts)
require(ggplot2)
mongo=mongo.create()
 
buf = mongo.bson.buffer.create()
mongo.bson.buffer.append(buf, "tickerId", 20L)
mongo.bson.buffer.start.object(buf, "size")
mongo.bson.buffer.append(buf, "$exists", "true")
mongo.bson.buffer.finish.object(buf)
 
query = mongo.bson.from.buffer(buf)
 
count = mongo.count(mongo,'quotes.trinti',query)
cursor=mongo.find(mongo,'quotes.trinti',query)
 
#############  very slow code #############
#size=''
#system.time(
#while(mongo.cursor.next(cursor)){
#  temp=(mongo.cursor.value(cursor));
#  if(is.xts(size))
#    size=rbind(size,xts(cbind(mongo.bson.value(temp,"field"),mongo.bson.value(temp,"size")),order.by=as.POSIXct(mongo.bson.value(temp,"tstamp")/1000,origin='1970-01-01',tz='Europa/Paris')))
#  else
#    size=xts(cbind(mongo.bson.value(temp,"field"),mongo.bson.value(temp,"size")),order.by=as.POSIXct(mongo.bson.value(temp,"tstamp")/1000,origin='1970-01-01',tz='Europa/Paris'))
})
#############  end very slow  #############
 
size=matrix(nrow=count,ncol=3)
counter=1
system.time(
  while(mongo.cursor.next(cursor))
  {
    temp=(mongo.cursor.value(cursor));
    size[counter,1]=mongo.bson.value(temp,"field");
    size[counter,2]=mongo.bson.value(temp,"size");
    size[counter,3]=mongo.bson.value(temp,"tstamp");
    counter=counter+1;
    if(counter>count)break;
    })
size=xts(size[,1:2],order.by=as.POSIXct(size[,3]/1000,origin='1970-01-01',tz='Europe/Paris'))
colnames(size)=c('field','size')
 
 
buf = mongo.bson.buffer.create()
mongo.bson.buffer.append(buf, "tickerId", 26L)
mongo.bson.buffer.start.object(buf, "price")
mongo.bson.buffer.append(buf, "$exists", "true")
mongo.bson.buffer.finish.object(buf)
 
query = mongo.bson.from.buffer(buf)
count = mongo.count(mongo,'quotes.trinti',query)
 
cursor=mongo.find(mongo,'quotes.trinti',query)
price=matrix(nrow=count,ncol=3)
counter=1
system.time(
  while(mongo.cursor.next(cursor))
  {
    temp=(mongo.cursor.value(cursor));
    price[counter,1]=mongo.bson.value(temp,"field");
    price[counter,2]=mongo.bson.value(temp,"price");
    price[counter,3]=mongo.bson.value(temp,"tstamp");
    counter=counter+1;
    if(counter>count)break;
  })
price=xts(price[,1:2],order.by=as.POSIXct(price[,3]/1000,origin='1970-01-01',tz='Europe/Paris'))
price=(price[which(price[,2]>0)])
 
colnames(price)=c('field','price')
 
quotes=cbind(price[,2][price[,1]==1],
             #cac40.volume[,2][cac40.volume[,1]==0],
             price[,2][price[,1]==2],
             #cac40.volume[,2][cac40.volume[,1]==3],
             price[,2][price[,1]==4]
             ,size[,2][size[,1]==5]
             )
 
quotes[,1]=na.locf(quotes[,1])
quotes[,2]=na.locf(quotes[,2])
quotes[,3]=na.locf(quotes[,3])
quotes[which(is.na(quotes[,4])),3]=NA
 
temp=tail(head(quotes,3000),1000)
temp=data.frame(ind=1:NROW(temp),trd=as.numeric(temp[,3])                
                ,bid=as.numeric(temp[,1]),ask=as.numeric(temp[,2])
                ,size=as.numeric(temp[,4])
                )
temp=melt(temp,id=c('ind'),na.rm=TRUE)
x=temp[which(temp$variable=='trd'),]
 
rez=temp[which(temp$variable!='trd'),]
rez=rez[which(rez$variable!='size'),]
a=temp[which(temp$variable=='size'),][,3]
ggplot(rez,aes(x=ind,y=value,color=variable))+geom_line()+geom_point(data=x,aes(size=a))

Today is Lithuania’s independence day and I have created R user group in Lithuania – VilniusR. If you are near by please follow the link, sign up and I hope that we will have a meeting soon.

In my previous post, I tried to show, that Rcpp is 1000 faster than pure R and that generated the fuss in the comments. Being lazy, I didn’t vectorize R code and at the end I was comparing apples vs oranges.

To fix that problem, I built a new script, where I’m trying to compare apples against apples. First piece of code named “ifelse R” uses R “ifelse” function to vectorize code. Second piece of code is fully vectorized code written in R, third – pure C++ code and the last one is C++, where  Rcpp ”ifelse” function is used.

Here we go – vectorization truly helps, but pure C++ code still 23 times faster. Of course you pay the price when writing it in C++.
I found a bit strange, that vectorized C++ code doesn’t perform that well…

You can get the code from github or review it below:

#Author Dzidorius Martinaitis
#Date 2012-02-01
#Description http://www.investuotojas.eu/2012/02/01/vectorized-r-vs-rcpp
 
bid = runif(50000000,5,9)
ask = runif(50000000,5,9)
close = runif(50000000,5,9)
 
x=data.frame(bid=bid,ask=ask,last_price=close)
rez=0
 
###########    ifelse R  #################
answ=as.vector(system.time(
{
rez = ifelse(x$last_price>0,ifelse(x[, "bid"] > x[, "last_price"], x[, "bid"], ifelse((x[, "ask"] > 0) & (x[, "ask"] < x[, "last_price"]), x[, "ask"], x[, "last_price"])), 0.5*(x[, "ask"] + x[,"bid"]))
})[1])
###########   end ifelse R  #################
 
###########    vectorized R  #################
 
answ=append(answ,system.time(
{
lgt0 = x$last_price > 0
bgtl = x$bid > x$last_price
agt0 = x$ask > 0
altl = x$ask > x$last_price
rez = x$last_price
rez[lgt0 & agt0 & altl] = x$ask[lgt0 & agt0 & altl]
rez[lgt0 & bgtl] = x$bid[lgt0 & bgtl]
rez[!lgt0] = (x$ask[!lgt0]+x$bid[!lgt0])/2
}
)[1])
###########   end vectorized R  #################
 
#C++ code starts here
 
library(inline)
library(Rcpp)
 
###########    pure C++  #################
 
code='
NumericVector bid(bid_);NumericVector ask(ask_);NumericVector close(close_);
int bid_size = bid.size();
NumericVector ret(bid_size);
for(int i =0;i<bid_size;i++)
{
  if(close[i]>0)
  {
    if(bid[i]>close[i])
    {
      ret[i] = bid[i]; 
    }
    else if(ask[i]>0 && ask[i]<close[i])
    {
      ret[i] = ask[i];//
    }
    else
    {
      ret[i] = close[i];//
    }
  }
  else
  {
    ret[i]=(bid[i]+ask[i])/2;
  }
 
}
return ret;
'
getLastPrice <- cxxfunction(signature( bid_ = "numeric",ask_ = "numeric",close_="numeric"),body=code,plugin="Rcpp")
rez=0
answ=append(answ,system.time(
  {
    rez=getLastPrice(as.numeric(x$bid),as.numeric(x$ask),as.numeric(x$last_price))
  })[1])
 
###########   end pure C++  #################
 
#summary(rez)
 
 
###########    vectorized C++  #################
code='
NumericVector bid(bid_);NumericVector ask(ask_);NumericVector close(close_);
int bid_size = bid.size();
NumericVector ret=ifelse(close>0,ifelse(bid >close, bid, ifelse(ask > 0,ifelse(ask < close,ask, close),close)), 0.5*(ask + bid));
return ret;
'
getLastPrice <- cxxfunction(signature( bid_ = "numeric",ask_ = "numeric",close_="numeric"),body=code,plugin="Rcpp")
rez=0
answ=append(answ,system.time(
{
  rez=getLastPrice(as.numeric(x$bid),as.numeric(x$ask),as.numeric(x$last_price))
}
)[1])
 
###########   end vectorized C++  #################
 
#summary(rez)
names(answ)=c('ifelse R','vectorized R','pure C++','vectorized C++')
 
library(ggplot2)
a=data.frame(ind=1:4,val=answ)
ggplot(a,aes(ind,val))+geom_point(legend=F)+geom_text(aes(label=names(answ),hjust=c(-0.2,-0.2,-0.2,0.8),vjust=c(0,0,0,-1)),size=4)

While ago I built two R scripts to track OMX Baltic Benchmark Fund against the index. One script returns the deviation of  fund from the index and it works fast enough. The second calculates the value of the fund every minute and it used to take for while. For example, it spent 2 minutes or more to get the values for one day. Here is an example of the result:

Following piece of code was in question:

for(y in 1:NROW(x))
 {
    z=x[y,]
    if(as.numeric(z$last_price>0))
    {
      if(as.numeric(z$bid>z$last_price))rez[y]=z$bid
      else if(as.numeric(z$ask)>0 &amp; as.numeric(z$ask)<z$last_price)rez[y]=z$ask
      else rez[y]=z$last_price
    }
    else
    {
      rez[y]=(z$ask+z$bid)/2
    }
 }

The code above loops over time series and based on set of rules tries to decide which price (bid, ask or previous one) to use for calculations. Pure R script used to take 100 seconds to derive the price.

During the weekend I found time to watch very interesting Rcpp presentation. To my surprise, there are numerous ways to seamlessly integrate C++ into R code. So, I decided to rewrite the code above in C++ (Rcpp and inline packages were used).

#c++ code embed in code value
code='
NumericVector bid(bid_);NumericVector ask(ask_);NumericVector close(close_);NumericVector ret(ask_);
int bid_size = bid.size();
for(int i =0;i<bid_size;i++)
{
  if(close[i]>0)
  {
    if(bid[i]>close[i])
    {
      ret[i] = bid[i];
    }
    else if(ask[i]>0 &amp;&amp; ask[i]<close[i])
    {
      ret[i] = ask[i];//
    }
    else
    {
      ret[i] = close[i];//
    }
  }
  else
  {
    ret[i]=(bid[i]+ask[i])/2;
  }
 
}
return ret;
'
#a glue function between C++ and R
getLastPrice = cxxfunction(signature( bid_ = "numeric",ask_ = "numeric",close_="numeric"),body=code,plugin="Rcpp")
 
#and the call of the function
getLastPrice(as.numeric(x$bid),as.numeric(x$ask),as.numeric(x$last_price))

What did I get in return? Well, 0.1 of a second instead of 100 seconds!

During the summer I was contacted by a hedge fund from Bahamas. The fund was looking for someone with R language skills on-site and insisted for phone interview. Besides obvious questions about finance, statistics, coding and how many tennis balls can fit in Boeing 747 (ok, this question was omitted), they wanted to know if I code in C++. So, I told them true – the last time I wrote a line in C++ was more or less 10 years ago. Long story short – it made me thinking about existence of C++.

10 yeas ago I was told, that C++ is going to disappear soon and Java is the king. At that time neither HN nor stackoverflow existed (meaning, that I had to rely on limited source), so I took it for granted, so here I am.

What do we have 10 years later? Neither C++ is dead, nor Java is sexy anymore. Actually is opposite – if you use Java, then you are clumsy programmer with lack of imagination. Does it sounds offensive? Then read for example the comments of Scala vs Java article and you will get the same feeling. Replacement of  Sun with Oracle does not help either.

But lets go back to C++. Google trends says, that C++ enjoys either maturity or decline. However, if you concentrate on specific industry, you will have a different picture. Kernel development (C not C++), game industry, number crunching, data mining, finance – where C++ matters. I know, I know, that you can write a magic code with Ruby or Python and it will perform almost as C++. And I saw a video, where guys were claiming, that they tuned “a bit”  Java and now it is able to deal with more that 1 million requests a second. Only the thing they did was elimination of garbage collection. The question – is it really worth of doing that way?

Next thing is to check what is demand for C++. Time to time I scroll through HN to be millionaires list and strangely enough C++ is demanded for back end systems, where performance or data amount is an issue. However, if you are targeting finance industry exclusively, then you may find this discussion interesting. Basically, it says, that there is stable demand for C++. Worth to say nonetheless, that C++ is mostly used by front and middle offices (where quants live) and the demand diminish in back office.

With such subjective study in mind I purchased C++ Primer by Stanley B. Lippman, which I recommend to beginners or disillusioned users like me. So far I built 2 small projects and in which one of them I parse 1.5 million terms to get a list of most used terms. R language does that in minutes, C++ in seconds.

Welcome back C++. It seems, that I miss you.

When dealing with transaction cost analysis, a stock’s volume is assumed to be stable or foreseeable.  However, there is different picture, then we are dealing with an illiquid stock.

It is relatively easy to forecast the volume of a liquid stock, because trading volume has high autocorrelation – the volumes at t and t+1 are correlated. For example, let’s take a look at BPN Paribas stock  autocorrelation figure:

X axis shows the lag between the days and Y axis shows percentage of the correlation. For BNP Paribas stock we have 60 % correlation between t and t+1 and 30 % between t and t+2.

Now, let’s look what is autocorrelation  of an illiquid stock:

The figure above shows, that autocorrelation for this emerging market stock is zero. That means, that we can’t forecast tomorrow’s volume, based on today’s volume.

Imagine, that a portfolio manager has to liquidate 100 000 of illiquid stock where the median of daily volume is 4000 and participation rate has to be maximum 20 %. How many days it will take to liquidate the position?

Because the daily volume of the stock is very volatile, we need more randomness in our forecast. To do that, we can use bootstrap – let’s take the last 250 data-points of the volume and generate 10 000 or 100 000 new time series. Once the have a bunch of new time series, let’s check how many days it would take to liquidate 100 000 stock position in each. Finally, we collect the numbers of the days needed for liquidation  and form a new vector. The result is following:

The histogram above shows, that it will take maximum 80 days to liquidate 100 000 of illiquid stock with 95% confidence.

#volume - vector of the stock's volume
#shareNumber - number of the shares to liquidate
#loop - number of the time-series to be created
#participationRate - what is going to be participation rate 
simNumberOfDays<-function(volume,shareNumber,loop=10000,participationRate=0.05)
{
  rez=matrix(nrow=loop)
  for(i in 1:loop)
  {
    x=nasa[sample(NROW(volume))]
    y=cumsum(x*participationRate)
    rez[i,]=which(y>shareNumber)[1]
  }
  return(rez)
}
rez=simNumberOfDays(illiquidStock,100000,10000,0.2)