#Simulation of Gaussian stochastic processes and Gaussian random fields #as high dimensional Gaussian random vectors library(RandomFields) #Plot of covariance functions of the Whittle-Matern class x <- 0.01*(0:500) plot(x,exp(-(x/2)^2),xlab="dist",ylab="covariance",type="l") a <- c(0.25,0.5,0.75,1,2,20) for (i in 1:length(a)){ cv <- CovarianceFct(x=x,model="whittle", param=c(NA,1,0,1/sqrt(a[i]),a[i])) lines(x,cv,col=i+1) } abline(h=0) abline(v=0) # Main difference is the behavior for small lags, decay to zero is similar #Simulation of such processes #The method "direct" uses the Cholesky decomposition #The method "circulant" uses circular embedding and is faster par(mfrow=c(3,1)) a <- 0.5 x <- GaussRF(x=c(0,20,0.01),grid=TRUE,gridtriple=TRUE, model="whittle",param=c(0,1,0,1/sqrt(a),a),method="circulant") plot(0.01*(1:1000),x[1:1000],type="l",ylim=range(x),xlab="time",ylab="x") plot(0.01*(1001:2000),x[1001:2000],type="l",ylim=range(x),xlab="time",ylab="x") acf(x,lag.max=500,demean=F) rho <- CovarianceFct(x=0.01*(0:500),model="whittle", param=c(NA,1,0,1/sqrt(a),a)) lines(0:500,rho,col=2) #Changing of a changes the local smoothness #Simulation of random fields (= random functions of 2 or more variables) #with these covariance functions. par(mfrow=c(2,2)) for (a in c(0.25,0.5,1,2)) { x <- GaussRF(seq(0,10,0.05),seq(0,10,0.05),grid=T,model="whittle", param=c(0,1,0,1/sqrt(a),a),method="circulant") image(x,col = terrain.colors(100)) } #Brownian motion and Brownian Bridge #Simulation as a random walk (independent increments) f.sim1.BM <- function(T=1,length=1000,start=0) { ## Purpose: Simulation of Brownian motion on a regular lattice of [0,T] ## ---------------------------------------------------------------------- ## Arguments: length = number of time points - 1 ## ---------------------------------------------------------------------- ## Author: Hans-Ruedi Kuensch, Date: 1 Jun 2006, 16:16 c(start, start + sqrt(T/length)*cumsum(rnorm(length))) } par(mfrow=c(2,1)) plot(seq(0,1,length=1001),f.sim1.BM(length=1000), type="l",ylim=c(-3,3),xlab="t",ylab="B(t)") for (i in (2:5)) lines(seq(0,1,length=1001),f.sim1.BM(length=1000),col=i) #Simulation according to the Lemma in the script: f.sim0.BM <- function(n=8) { ## Purpose: Simulation of Brownian motion at time poitns j*2^{-n} T, ## j=0,1,2,..., 2^n ## ---------------------------------------------------------------------- ## Arguments: ## ---------------------------------------------------------------------- ## Author: Hans-Ruedi Kuensch, Date: 1 Jun 2006, 16:16 k <- 2^n y <- matrix(0,nrow=k+1,ncol=n) y[,1] <- rnorm(1)*(0:k)/k for (j in (1:(n-1))) { x <- rnorm(2^(j-1)) int <- (1:2^(n-j))/2^(n-j) z <- 0 for (i in (1:length(x))) z <- c(z,x[i]*int,x[i]*(1-int)) y[,j+1] <- z/sqrt(2^(j+1)) } y } par(mfrow=c(3,3)) n <- 9 time <- (0:(2^n))/2^n y <- f.sim0.BM(n=n) B <- apply(y,1,sum) for (i in (1:n)) { plot(time,y[,i],ylim=range(c(B,y)),type="l",ylab="y",col=4) if (i>1) lines(time,apply(y[,1:i],1,sum)) } par(mfrow=c(2,1)) n <- 10 time <- (0:(2^n))/2^n B <- apply(f.sim0.BM(n=n),1,sum) plot(time,B,type="l",ylim=c(-2.5,2.5)) for (i in (2:5)) { B <- apply(f.sim0.BM(n=n),1,sum) lines(time,B,col=i) } # Simulation using eigenvalues and eigenvectors # (called Karhunen-Loeve-expansion) psi <- function(i,t) sqrt(2)*sin((i+0.5)*pi*t) #Eigenvectors time <- seq(0,1,length=1001) plot(time,psi(0,time),type="l",ylim=c(-1.5,1.5),ylab="Eigenfunctions") for (i in (1:5)) lines(time,psi(i,time),col=i+1) plot(time,psi(6,time),type="l",ylim=c(-1.5,1.5),ylab="Eigenfunctions") for (i in (1:5)) lines(time,psi(i+6,time),col=i+1) n <- 10000 # Number of eigenfunctions used lambda <- 1/(pi*((0:n)+0.5))^2 # square root of eigenvalues plot(1:n,lambda[-1],type="h",log="xy") f.sim2.BM <- function(z,time=seq(0,1,length=1001)) { ## Purpose: Simulation of Brownian motion using Karhunen-Loeve expansion ## ---------------------------------------------------------------------- ## Arguments: z=Koeffizienten in der Entwicklung (i.i.d N(0,1)) ## ---------------------------------------------------------------------- ## Author: Hans-Ruedi Kuensch, Date: 1 Jun 2006, 16:16 n <- length(z) -1 lambda <- 1/(pi*((0:n)+0.5)) drop((z*lambda)%*%outer(0:n,time,FUN="psi")) } n <- 1000 z <- rnorm(n+1) plot(time,f.sim2.BM(z=z,time=time),type="l",ylab="B") lines(time,f.sim2.BM(z=z[1:11],time=time),col=2) lines(time,f.sim2.BM(z=z[1:51],time=time),col=4) lines(time,f.sim2.BM(z=z[1:101],time=time),col=6) #Brownian Bridge (not done in the course) f.sim.BB <- function(start=0,end=0,time=seq(0,1,by=0.001)) { ## Purpose: Simulation einer Brown'schen Bruecke. ## ---------------------------------------------------------------------- ## Arguments: time soll ein geordneter Vektor der Laenge >=3 sein. ## start und end sind die Werte der Brown'schen Bewegung zu ersten ## und zur letzten Zeit, die Funktion simuliert an den Zwischenzeiten. ## ---------------------------------------------------------------------- ## Author: Hans-Ruedi Kuensch, Date: 2 Jun 2006, 10:30 n <- length(time) z <- rnorm(n-2) d.t <- time[2:(n-1)]-time[1:(n-2)] l.t <- time[n]-time[1:(n-2)] u.t <- d.t/l.t v.t <- 1-u.t z <- u.t*end + sqrt(d.t*v.t)*z b <- c(rep(start,n-1),end) for (i in (2:(n-1))) b[i] <- v.t[i-1]*b[i-1]+z[i-1] b} plot(time,f.sim.BB(time=time),type="l",ylab="B") plot(time,f.sim.BB(end=rnorm(1),time=time),type="l",ylim=c(-3,3),ylab="B") for (i in (2:5)) lines (time,f.sim.BB(end=rnorm(1),time=time),col=i) abline(h=0)