######## 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")