josef teichmann :: personal homepage at ETH Zurich

Probability and Statistics FS 2010

We mainly use R for numerical implementations, code can be found here. You can install R on your computer via CRAN. Load the code into R and run the function with the appropriate parameters fixed. Usually default values of parameters are provided.

Numerical and MC-integration of a 1d integral

Simulation of numbers of long runs for random walks and their statistics

Simulation and (beautiful) visualization of stochastic 3d-population dynamics according to a Wishart process

# representation as heat map or surface (uncomment it and comment the heat map!) of one trajectory of the time evolution of a 3d continuous matrix valued affine process 
matrix.valued.affine.process3d <- function(X0=diag(3,3),b=diag(c(0,0,0)),beta=diag(3,3),alpha=diag(3,3),n=50,time=1,m=50)
{
y<-seq(0,2*pi,length.out=(2*m)) # initialization of two vectors representing the angles on a 2d sphere
x<-seq(0,pi,length.out=m)
z<-matrix(nrow=m,ncol=(2*m)) # initialization of a matrix corresponding to x[i] and y[j]
for (i in (1:n)) # outer loop along for the euler scheme
{
X<-X0
Y<-eigen(X)
Z<-Y$vectors
sqrtX<-Z %*% sqrt(abs(diag(Y$values))) %*% t(Z)
W<-matrix(3,3,data=rnorm(9))
X<-X+(b+beta %*% X+X %*% t(beta))*(time/n)+ sqrtX %*% (W *(sqrt(time/n))) %*% alpha + t(alpha) %*% (t(W) *(sqrt(time/n))) %*% sqrtX
for (i in (1:m)) # inner loop for the evaluation of the time i result of X on the 2d-sphere
for (j in (1:(2*m)))
{{
z[i,j]<-t(c(cos(i*pi/m)*cos(j*pi/m),cos(i*pi/m)*sin(j*pi/m),sin(i*pi/m))) %*% X %*% t(t(c(cos(i*pi/m)*cos(j*pi/m),cos(i*pi/m)*sin(j*pi/m),sin(i*pi/m)))) 
}}
#persp(x,y,z,theta = 30, phi = 30, expand = 0.5, col = "lightblue",box=FALSE) # drawing of a perspective plot in lightblue without box
filled.contour(x,y,z,axes=FALSE) # drawing of a heat map without axes
}
}

Determinant process for a 2d-Wishart process

# the function calculates m timeseries of a the determinant of a 2d-matrix-valued affine diffusion process

det.process.2d <- function(X0=diag(c(1,1)),b=diag(c(0.4,0.1)),beta=diag(c(-0.1,-0.2)),alpha=diag(c(0.3,0.1)),n=500,time=1,m=10)
{
x<-(1:657) # R-colors in numbers
y<-sample(x) # a random sample of x
for (j in (1:20)) # the loop for the m timeseries
{
X<-X0 # initialization of X variable
ts<-vector(length=n) # initialization of the ts variable
for (i in (1:n)) # inner loop along the euler discretization
{
Y<-eigen(X) # eigenvalues and vectors of X
Z<-Y$vectors
sqrtX<-Z %*% sqrt(diag(abs(Y$values))) %*% t(Z) #calculation of the squareroot of X
W<-matrix(2,2,data=rnorm(4)) # drawing of the random matrix (not symmetric) W with variance 1
X<-X+(b+beta %*% X+X %*% t(beta))*(time/n)+ sqrtX %*% (W *(sqrt(time/n))) %*% alpha + t(alpha) %*% (t(W) *(sqrt(time/n))) %*% sqrtX
#ts[i]=sum(diag(X))
ts[i]=det(X)
}
if (j==1) plot(ts,type="l",ylim=c(0,5)) # plot the first time series
else lines(ts,col=colors()[y[j]]) # add all additional ones with a randomly chosen color
}
}