#Week 2:
# Simulations from simple stochastic models
# New functions: filter

#Simulation of a (linear) moving average
coeff <- rep(1,5) #different choices of the parameters
coeff <- c(1,1.9,1)
coeff <- c(1,-2,1)
eps <- ts(rnorm(200))
x <- filter(eps,filter=coeff,sides=1)
plot(ts.union(eps,x))
help(filter) #explains the function filter


#Simulation of (linear) autoregressive processes
phi <- 0.5 #different choices of the parameter
phi <- -0.5
phi <- 2
phi <- 0.98
x <- filter(ts(rnorm(200)),filter=phi,method="recursive",init=rnorm(1))
#different initial value
x <- filter(ts(rnorm(200)),filter=phi,method="recursive",init=-10)
plot(x)
#same thing without the function filter
x <- ts(rnorm(200))
for (i in (2:100)) x[i] <- phi*x[i-1]+x[i]

#Simulation of nonlinear autoregressions
f.1 <- function(x) x*(0.5+25/(1+x^2)) #Example of Andrade-Netto
n <- 1000
sigma <- sqrt(10) 
x <- ts(rnorm(n))
for (i in (2:n)) x[i] <- f.1(x[i-1]) + sigma*x[i]
plot(x)

z <- seq(-20,20,by=0.01) #plotting the regression function
plot(z,f.1(z),type="l")
abline(a=0,b=1)
points(x[-1000],x[-1])

f.2 <- function(x) 3.8*x*(1-x) #logistic growth
n <- 500
sigma <- 0.01 # s.d. of the noise. See what happens if sigma=0
x <- ts(rnorm(500))
x[1] <- 0.5
for (i in (2:500)) {
    x[i] <- f.2(x[i-1]) + sigma*x[i]
    x[i] <- max(0,min(1,x[i])) #keep the process in [0,1]
  }
plot(x)

z <- seq(0,1,by=0.01)
plot(z,f.2(z),type="l")
abline(a=0,b=1)
points(x[-n],x[-1])
hist(x,nclass=20)


# Lorenz model (3-dim. differential equation plus noise)
f.lorenz <- 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.05 #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*f.lorenz(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")

#ARCH-processes
n <- 1000
phi <- 0.9
eps <- ts(rnorm(n))
x <- eps
for (i in (2:n)) x[i] <- sqrt(1+phi*x[i-1]^2)*eps[i]
plot(ts.union(eps,x),plot.type="single",col=2:1)
abline(h=0)