ESSA2013 Conference

It has been just announced that during ESSA2013 conference I am planning to organize a special track on "Statistical analysis of simulation models". I hope to get some presentations using GNU R to promote it in social simulation community.

It is obvious that GNU R excels in analysis of simulation data. However, very often it can be neatly used to implement simulations themselves.

For instance I have recently implemented a simulation model proposed in Section 4 of Volatility Clustering in Financial Markets: Empirical Facts and Agent–Based Models paper by Rama Cont. The model is formulated as follows (I give only its brief description; please refer to the paper for more details).

Consider market with n trading agents and one asset. We simulate market for times periods. In each period each agent can buy asset, sell it or do nothing.
Asset return r[i] in period i equals to number of buy orders minus number of sell orders divided by number of agents n and multiplied by normalizing constant max.r. Thus it will always lie in the interval [-max.r,max.r].
Agents make buy and sell decisions based on random public information about an asset. The stream of signals are IID normal random variables with mean 0 and standard deviation signal.sd. Each investor holds an internal non negative decision making threshold. If signal is higher than threshold level buy decision is made. If it is lower than minus threshold level asset is sold. If signal is not strong enough investor does nothing.
After return r[i] is determined each investor with probability p.update performs threshold update to abs(r[i]).
As you can see the description is quite lengthily. However, the implementation of the model in GNU R is a genuine snippet as can be seen below:

cont <- function(times, n, signal.sd, max.r ,p.update) {

      threshold <- vector("numeric", n)

      signal <- rnorm(times, 0, signal.sd)

      r <- vector("numeric", times)

      for (i in 1:times) {

            r[i] <- max.r * (sum(signal[i] > threshold)

                         - sum(signal[i] < (-threshold))) / n

            threshold[runif(n) < p.update] <- abs(r[i])

      }

      return(r)

}

And an additional benefit is that one can analyze the simulation results in GNU R also. Here is a very simple example showing the relationship between signal.sd and standard deviation of simulated returns (the initial burn in period in the simulation is discarded):

cont.sd <- function(signal.sd) {

    sd(cont(10000, 1000, signal.sd, 0.1, 0.05)[1000:10000])

}

sd.in <- runif(100, 0.01, 0.1)

sd.out <- sapply(sd.in, cont.sd)

plot(sd.in,sd.out)

and here is the resulting plot: