Prof Darren Wilkinson Professor of Stochastic Modelling
School of Mathematics & Statistics
Newcastle University

Markov chain simulation and analysis

This week we will look at how to simulate Markov chains, and in particular, how to gain an empirical understanding of their stationary distributions. First read carefully through the following example, trying it out as you go along, then tackle the exercises below.

Exercises

1.
Try simulating AR(1) processes with different parameter values. Note that you will have problems with parameter values outside the range [-1,1]! Try starting off with parameters inside this range, and see what happens as you get close to the boundaries. eg. try 0.9, 0.95, 0.99 on the positive side, and -0.9, -0.95, -0.99 on the negative side. A parameter value of -1 often produces an interesting plot! Notice how the chain converges much more quickly for some parameter values than others. Note also how the location, scale and shape of the stationary distribution changes as the parameter values change.
2.
What is the theoretical asymptotic mean and variance of the process, in terms of the parameter, a? Make sure you have the right answer by comparing with your empirical findings for various values of a.
3.
What is the (approximate) standard error of the sample mean as a function of n and a? Hint: remember that the observations are not independent! Use the information in your lecture notes.
4.
For a particular sample of 10,000 that you generate for a=0.995, construct a confidence interval for the mean. Pretend that you don't know a, and use the sample lag 1 autocorrelation to estimate it. Compare that to the confidence interval you would have if the observations were independent.
5.
Calculate the sample lag 2 autocorrelation in the chain from the previous question, and use this to verify that the assumption of an underlying AR(1) process is reasonable.
6. (for postgrads)
Try out the C code for chain generation. See if you can use the GSL for variate generation. Try out the simple method for calling executables from R. Then try compiling your C code as a shared object, and dynamically loading it into R, and calling it directly.

Darren J Wilkinson Book: Stochastic Modelling for
	      Systems Biology  
darren.wilkinson@ncl.ac.uk
http://www.staff.ncl.ac.uk/d.j.wilkinson/