## Simulate Gaussian processes
# Set-up
L = 10 # upper limit for time say

## Exponential
set.seed(1134)
# Parameters
pE = sqrt(8*0.5)/3
sig = 1

# Grid
xx = seq(0, L, length.out = 1000)
h = abs(as.matrix(dist(xx)))

# Description
mu = rep(0, length(xx))
Sig = sig^2*exp(-pE*h)

# Simulate several realizations
nR = 3
yMat = matrix(0, nrow = nR, ncol = length(xx))
Lt = t(chol(Sig))
for(i in 1:nR){
  yMat[i,] = Lt%*%rnorm(length(xx))
}
yAll = yMat[1,]

# Plot
matplot(xx, t(yMat), col=1:nR, type="l", lty=1,
        xlim = c(xx[1], tail(xx,1)), ylim = c(-4*sig, 4*sig), 
        xlab = "Time", ylab = "Response", cex.axis = 1.5, 
        cex.lab = 1.5, lwd = 1.5,
        main = "Exponential covariance function")


## Matern-like
set.seed(1134)
# Parameters
pM = sqrt(8*1.5)/3
sig = 1


# Description
mu = rep(0, length(xx))

Sig = sig^2*(1+pM*h)*exp(-pM*h)

# Simulate several realizations
nR = 3
yMat = matrix(0, nrow = nR, ncol = length(xx))
Lt = t(chol(Sig))
for(i in 1:nR){
  yMat[i,] = Lt%*%rnorm(length(xx))
}
yAll = rbind(yAll, yMat[1,])

# Plot
matplot(xx, t(yMat), col=1:nR, type="l", lty=1,
        xlim = c(xx[1], tail(xx,1)), ylim = c(-4*sig, 4*sig), 
        xlab = "Time", ylab = "Response", cex.axis = 1.5, 
        cex.lab = 1.5, lwd = 1.5,
        main = "Matérn-type covariance function")


## Gaussian
set.seed(1134)
# Parameters
pG = 0.22
sig = 1

# Description
mu = rep(0, length(xx))
Sig = sig^2*exp(-h^2*pG)

# Simulate several realizations
nR = 3
yMat = matrix(0, nrow = nR, ncol = length(xx))
Lt = t(chol(Sig+diag(length(xx))*1e-10))
for(i in 1:nR){
  yMat[i,] = Lt%*%rnorm(length(xx))
}
yAll = rbind(yAll, yMat[1,])

# Plot
matplot(xx, t(yMat), col=1:nR, type="l", lty=1,
     xlim = c(xx[1], tail(xx,1)), ylim = c(-4*sig, 4*sig), 
     xlab = "Time", ylab = "Response", cex.axis = 1.5, 
     cex.lab = 1.5, lwd = 1.5,
     main = "Gaussian covariance function")

## One of realisation for each covariance function
matplot(xx, t(yAll), 
        col=1:3, lwd=1.5, type="l", lty=1, 
        xlim = c(xx[1], tail(xx,1)), 
        ylim = c(-4*sig, 4*sig), 
        xlab = "Time", ylab = "Response", 
        cex.axis = 1.5, cex.lab = 1.5)

legend('topright', legend = c("Exponential", "Matérn-like", "Gaussian"), 
       lwd = 2, col = c("black", "red", "green"))

## Correlation functions
plot(xx, exp(-pE*abs(xx)), xlab = "Time difference", ylab ="Correlation", lwd = 2, type = "l")
lines(xx, (1+pM*xx)*exp(-pM*xx), lwd = 2, col = "blue", lty = 2)
lines(xx, exp(-pG*xx^2), lwd = 2, col = "red", lty = 3)
#abline(h = 0, lwd = 2)   
legend('topright', legend = c("Exponential", "Matérn-like", "Gaussian"), lwd = 2, col = c("black", "blue", "red"), 
       lty = c(1,2,3))