######## GAUSSIAN WHITE NOISE N <- 300 W <- ts(rnorm(N)) plot(W,type="b",xlab="time") plot(X,xlab="time",type="n",col=rgb(0.2,0.2,0.2),ylim=c(-50,50)) ######## RANDOM WALK for (sim in 1:100){ W <- rnorm(1000) X <- cumsum(W) lines(X,xlab="time",type="l",col=sim) } plot(exp(X/100),type="l") ## equivalently X[1] <- W[1] for (i in 2:N){ X[i] <- X[i-1] + W[i] } plot(X,xlab="time",type="b") ## equivalently X <- filter(W, filter=1, method="recursive") plot(X,xlab="time",type="b") ######### #(linear) moving average coeff <- c(1,-1) #different choices of the parameters coeff <- rep(1,5) X <- filter(W,filter=coeff,sides=1) plot(ts.union(W,X),type="b") plot(X[1:300],type="b", xlab="time") help(filter) #explains the function filter ##equivalently Z <- numeric(N) Z[1] <- W[1] for (i in 2:N){ Z[i] <- sum(W[i:max(1,i-4)]) } lines(Z,col=2) ######## linear auto-regressive coeff <- c(0.4, 0.2, -0.3) X <- filter(W,filter=coeff,method="recursive") plot(ts.union(W,X),type="b") plot(X,type="l") coeff <- c(1.99,-0.991) X <- filter(W,filter=coeff,method="recursive") plot(ts.union(W,X),type="b") plot(X,type="l") coeff <- c(-1.99,-0.991) X <- filter(W,filter=coeff,method="recursive") plot(ts.union(W,X),type="b") plot(X,type="l") ## without using the function filter X <- numeric(N) X[1] <- W[1] phi <- 0.99 for (i in (2:N)) X[i] <- phi*X[i-1]+W[i] plot(X,xlab="time",type="b") ### ARCH-process N <- 2000 phi <- 0.95 W <- ts(rnorm(N)) X <- numeric(N) X[1] <- W[1] for (i in (2:N)) X[i] <- sqrt(1+phi*X[i-1]^2)*W[i] plot(ts.union(W,X)) plot(X,type="l",xlab="time") plot(exp(cumsum(X)/100),type="l") ### Simulation of nonlinear autoregressions f <- function(x) x*(0.5+25/(1+x^2)) #Example of Andrade-Netto N <- 1000 sigma <- sqrt(10) W <- ts(rnorm(N)) X <- numeric(N) X[1] <- 0 for (i in (2:N)) X[i] <- f(X[i-1]) + sigma*W[i] plot(X) z <- seq(-20,20,by=0.01) #plotting the regression function plot(z,f(z),type="l",lwd=3) abline(a=0,b=1,col=2,lty=3,lwd=2) points(X[-N],X[-1]) f2 <- function(x) 3.8*x*(1-x) #logistic growth N <- 200 sigma <- 0.01 # s.d. of the noise. See what happens if sigma=0 W <- ts(rnorm(N)) X <- numeric(N) X[1] <- 0.5 for (i in (2:N)) { X[i] <- f2(X[i-1]) + sigma*rnorm(1) X[i] <- max(0,min(1,X[i])) #keep the process in [0,1] } plot(X,type="b") z <- seq(0,1,by=0.01) plot(z,f2(z),type="l",lwd=3) abline(a=0,b=1,col=2,lty=3) points(X[-N],X[-1]) hist(X,nclass=20) # Lorenz model (3-dim. differential equation plus noise) florenz <- function(x) { c(10*(x[2]-x[1]),28*x[1]-x[2]-x[1]*x[3],x[1]*x[2]-2.666667*x[3]) } N <- 10000 sigma <- 0.01 #standard deviation for one Euler step X <- ts(matrix(0,nrow=N+1,ncol=3),start=0,frequency=200,end=50) X[1,] <- c(-6,-5,25) #fixed starting point for (i in (1:N)) X[i+1,] <- X[i,] + 0.005*florenz(X[i,])+ sigma*rnorm(3) plot(X) par(mfrow=c(1,2)) plot(X[,1],X[,2],type="l") plot(X[,1],X[,3],type="l")