Prof Darren Wilkinson Professor of Stochastic Modelling

School of Mathematics & Statistics

Newcastle University

Simulating an AR(1) with exponential innovations

Let us consider simulating the AR(1) process {X(t)} defined by X(t) = a*X(t-1) +Z(t) where a is a constant and the {Z(t)} are independently distributed Exp(1) random variables. We will start the process off at X(1)=0.

Implementation in R

By now you should have a good idea how to do this.

arsim<-function (n, alpha) 
{
        vec <- vector("numeric", n)
        vec[1] <- 0
        for (i in 2:n) {
                vec[i] <- alpha * vec[i - 1] + rexp(1, 1)
        }
        vec
}

Note that since each observation depends on the previous one, there is no way to vectorise this procedure. Consequently, it is a little slow, and for simulating more complex Markov chains than this one, recoding into C once working may be desirable.

We can test it out and examine the results of the simulation as follows.

test<-arsim(10000,0.99)
par(mfrow=c(2,1))
plot(ts(test))
hist(test,30)

Here we create a vector test containing a simulated chain of length 10,000 with a parameter of 0.99 (quite close to 1). We then obtain a time series plot of the sequence, and a histogram of the observations.

You can see that the first few values should be discarded, as the chain takes a little while to settle down to its stationary distribution.

The R command cor(test[1:9999],test[2:10000]) gives a sample estimate of the lag 1 autocorrelation of the process. If the simulated chain is AR(1), then this is an estimate of the parameter of the process.

The R source code is available here: arsim.r

Implementation in C

An implementation in C is very straightforward. This code uses the GNU Scientific library (GSL) for random number generation.

int main(int argc, char *argv[])
{
  int i,n;
  float x,alpha;
  gsl_rng *r=gsl_rng_alloc(gsl_rng_mt19937); /* initialises GSL RNG */
  n=atoi(argv[1]);
  alpha=atof(argv[2]);
  x=0;
  for (i=0;i<n;i++)
    {
      x=alpha*x + gsl_ran_exponential(r,1);
      printf(" %2.4f \n",x);
    }
  return(0);
}

The full C source for this is here: arsim.c.

Now that we are dealing with procedures which cannot be vectorised in R, some people might be interested in how you can call C procedures from R as normal R functions. A simple way to do it is as follows. If your compiled C program is called a.out define the R function:

arsim_function(n,lambda)
{
        system(paste("a.out", n, lambda, ">tmp.dat"))
        scan("tmp.dat")
}

Now arsim can be called as before, but in fact the C program generates a data file which is then read by R.

In fact, this method is a bit clunky, unreliable and slow, so a better way to do it is to link your compiled C functions into the R executable. However, this is much easier to do under Unix than under Windows NT. There is some modified C code containing brief instructions how to do this under Unix here: arsim-link.c.

Other implementations

Those who are interested may also wish to look at corresponding code for LISP-STAT: arsim.lsp,Python: arsim.py and Sather: arsim.sa.

Darren J Wilkinson

darren.wilkinson@ncl.ac.uk

http://www.staff.ncl.ac.uk/d.j.wilkinson/

	Professor of Stochastic Modelling
	School of Mathematics & Statistics
	Newcastle University


		`darren.wilkinson@ncl.ac.uk`
		`http://www.staff.ncl.ac.uk/d.j.wilkinson/`