(This article was originally published at R snippets, and syndicated at StatsBlogs.)

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:

**Please comment on the article here:** **R snippets**