##################################################################################
# SPATIAL PROCESS: LAKES EXAMPLE IN LEGENDRE AND LEGENDRE (2012)
##################################################################################

###############################################
# 1) A very simple spatial process

s <- seq(1, 100, 1)
length(s)

Y <- rnorm(length(s), mean=0, sd=1)
plot(s, Y)

for (k in 1:1000)
	{
	Y <- rnorm(length(s), mean=0, sd=1)
	plot(s, Y, main = paste("Realization", k, " of the process"))
	Sys.sleep(1)
	}

#capture data generated by 1000 realizations of the process at selected spatial points 
y.20 <- rep(NA, 1000)
y.21 <- rep(NA, 1000)
y.40 <- rep(NA, 1000)
y.41 <- rep(NA, 1000)
y.60 <- rep(NA, 1000)
y.61 <- rep(NA, 1000)
y.80 <- rep(NA, 1000)
y.81 <- rep(NA, 1000)
 
for (k in 1:1000)
	{
	Y <- rnorm(length(s), mean=0, sd=1)
	plot(s, Y, main = paste("Realization", k, " of the process"))
	y.20[k] <- Y[20]
	y.21[k] <- Y[21]
	y.40[k] <- Y[40]
	y.41[k] <- Y[41]
	y.60[k] <- Y[60]
	y.61[k] <- Y[61]
	y.80[k] <- Y[80]
	y.81[k] <- Y[81]
	}

#examine the mean of the process at selected points (s.20, s.21, ..., s.81) 
selected.s <- list(y.20, y.21, y.40, y.41, y.60, y.61, y.80, y.81)
boxplot(selected.s, names=c("s.20", "s.21", "s.40", "s.41", "s.60", "s.61", "s.80", "s.81"), ylab="Y")

#examine the covariance of the process between selected pairs of spatial points
plot(y.20, y.21)
cov(y.20, y.21)
cor.test(y.20, y.21)

plot(y.40, y.41)
cov(y.40, y.41)
cor.test(y.40, y.41)

plot(y.60, y.61)
cov(y.60, y.61)
cor.test(y.60, y.61)

plot(y.80, y.81)
cov(y.80, y.81)
cor.test(y.80, y.81)

plot(y.20, y.40)
cov(y.20, y.40)
cor.test(y.20, y.40)

plot(y.20, y.60)
cov(y.20, y.60)
cor.test(y.20, y.60)

plot(y.20, y.80)
cov(y.20, y.80)
cor.test(y.20, y.80)

plot(y.40, y.60)
cov(y.40, y.60)
cor.test(y.40, y.60)

plot(y.40, y.80)
cov(y.40, y.80)
cor.test(y.40, y.80)

plot(y.60, y.80)
cov(y.60, y.80)
cor.test(y.60, y.80)


###############################################
# 2) A more complex spatial process: Y depends on X, and X depends on s

s <- seq(1, 100, 1)
length(s)

X <- s*0.1
plot(s, X)

YA <- X*0.5 + rnorm(length(s), mean=0, sd=1)
plot(X, YA)
plot(s, YA)

#true mean
E.YA <- mean(X*0.5)
E.YA

#sample mean
mean(YA)

#1000 realizations of the process
for (k in 1:1000)
	{
	YA <- X*0.5 + rnorm(length(s), mean=0, sd=1)
	plot(s, YA, main = paste("Realization", k, " of the process"))
	Sys.sleep(1)
	}

#capture data generated by 1000 realizations of the process at selected spatial points 
yA.20 <- rep(NA, 1000)
yA.21 <- rep(NA, 1000)
yA.40 <- rep(NA, 1000)
yA.41 <- rep(NA, 1000)
yA.60 <- rep(NA, 1000)
yA.61 <- rep(NA, 1000)
yA.80 <- rep(NA, 1000)
yA.81 <- rep(NA, 1000)
 
for (k in 1:1000)
	{
	YA <- X*0.5 + rnorm(length(s), mean=0, sd=1)
	plot(s, YA, main = paste("Realization", k, " of the process"))
	yA.20[k] <- YA[20]
	yA.21[k] <- YA[21]
	yA.40[k] <- YA[40]
	yA.41[k] <- YA[41]
	yA.60[k] <- YA[60]
	yA.61[k] <- YA[61]
	yA.80[k] <- YA[80]
	yA.81[k] <- YA[81]
	}

#examine the mean of the process at selected spatial points (s.20, s.21, ..., s.81) 
selected.s <- list(yA.20, yA.21, yA.40, yA.41, yA.60, yA.61, yA.80, yA.81)
boxplot(selected.s, names=c("s.20", "s.21", "s.40", "s.41", "s.60", "s.61", "s.80", "s.81"), ylab="Y")

#examine the covariance of the process between selected pairs of spatial points
plot(yA.20, yA.21)
cov(yA.20, yA.21)
cor.test(yA.20, yA.21)

plot(yA.40, yA.41)
cov(yA.40, yA.41)
cor.test(yA.40, yA.41)

plot(yA.60, yA.61)
cov(yA.60, yA.61)
cor.test(yA.60, yA.61)

plot(yA.80, yA.81)
cov(yA.80, yA.81)
cor.test(yA.80, yA.81)

plot(yA.20, yA.40)
cov(yA.20, yA.40)
cor.test(yA.20, yA.40)

plot(yA.20, yA.60)
cov(yA.20, yA.60)
cor.test(yA.20, yA.60)

plot(yA.20, yA.80)
cov(yA.20, yA.80)
cor.test(yA.20, yA.80)

plot(yA.40, yA.60)
cov(yA.40, yA.60)
cor.test(yA.40, yA.60)

plot(yA.40, yA.80)
cov(yA.40, yA.80)
cor.test(yA.40, yA.80)

plot(yA.60, yA.80)
cov(yA.60, yA.80)
cor.test(yA.60, yA.80)


################################################
# 3) Another spatial process: autocorrelation in X

s <- seq(1, 100, 1)
length(s)

E.YB <- E.YA

YB <- rep(NA, times=length(s))
YB

YB[1] <- E.YB + rnorm(1, mean=0, sd=1)
YB

YB[2] <- E.YB + (YB[1]-E.YB) + rnorm(1, mean=0, sd=1)
YB

for(i in 3:length(s))
	{
	YB[i] <- E.YB + (YB[i-1]-E.YB) + rnorm(1, mean=0, sd=1)
	}
YB
plot(s, YB)

#1000 realizations of the process
for (k in 1:1000)
	{
	YB <- rep(NA, times=length(s))
	#YB

	YB[1] <- E.YB + rnorm(1, mean=0, sd=1)
	#YB

	YB[2] <- E.YB + (YB[1]-E.YB) + rnorm(1, mean=0, sd=1)
	#YB

		for(i in 3:length(s))
		{
		YB[i] <- E.YB + (YB[i-1]-E.YB) + rnorm(1, mean=0, sd=1)
		}
	#YB
	plot(s, YB, main = paste("Realization", k, " of the process"))
	Sys.sleep(1)
	}


#capture data generated by 1000 realizations of the process at selected spatial points 
yB.20 <- rep(NA, 1000)
yB.21 <- rep(NA, 1000)
yB.40 <- rep(NA, 1000)
yB.41 <- rep(NA, 1000)
yB.60 <- rep(NA, 1000)
yB.61 <- rep(NA, 1000)
yB.80 <- rep(NA, 1000)
yB.81 <- rep(NA, 1000)
 
for (k in 1:1000)
	{
	YB <- rep(NA, times=length(s))
	#YB

	YB[1] <- E.YB + rnorm(1, mean=0, sd=1)
	#YB

	YB[2] <- E.YB + (YB[1]-E.YB) + rnorm(1, mean=0, sd=1)
	#YB

		for(i in 3:length(s))
		{
		YB[i] <- E.YB + (YB[i-1]-E.YB) + rnorm(1, mean=0, sd=1)
		}
	#YB
	plot(s, YB, main = paste("Realization", k, " of the process"))
	#Sys.sleep(1)

	yB.20[k] <- YB[20]
	yB.21[k] <- YB[21]
	yB.40[k] <- YB[40]
	yB.41[k] <- YB[41]
	yB.60[k] <- YB[60]
	yB.61[k] <- YB[61]
	yB.80[k] <- YB[80]
	yB.81[k] <- YB[81]
	}

#examine the mean of the process at selected spatial points (s.20, s.21, ..., s.81) 
selected.s <- list(yB.20, yB.21, yB.40, yB.41, yB.60, yB.61, yB.80, yB.81)
boxplot(selected.s, names=c("s.20", "s.21", "s.40", "s.41", "s.60", "s.61", "s.80", "s.81"), ylab="Y")

#examine the covariance of the process between selected pairs of spatial points
plot(yB.20, yB.21)
cov(yB.20, yB.21)
cor.test(yB.20, yB.21)

plot(yB.40, yB.41)
cov(yB.40, yB.41)
cor.test(yB.40, yB.41)

plot(yB.60, yB.61)
cov(yB.60, yB.61)
cor.test(yB.60, yB.61)

plot(yB.80, yB.81)
cov(yB.80, yB.81)
cor.test(yB.80, yB.81)

plot(yB.20, yB.40)
cov(yB.20, yB.40)
cor.test(yB.20, yB.40)

plot(yB.20, yB.60)
cov(yB.20, yB.60)
cor.test(yB.20, yB.60)

plot(yB.20, yB.80)
cov(yB.20, yB.80)
cor.test(yB.20, yB.80)

plot(yB.40, yB.60)
cov(yB.40, yB.60)
cor.test(yB.40, yB.60)

plot(yB.40, yB.80)
cov(yB.40, yB.80)
cor.test(yB.40, yB.80)

plot(yB.60, yB.80)
cov(yB.60, yB.80)
cor.test(yB.60, yB.80)


################################################
# 4) A spatial process in which Y depends on X, X depends on s, and there is autocorrelation in Y

s <- seq(1, 100, 1)
length(s)

X <- s*0.1
plot(s, X)

YC <- rep(NA, times=length(s))
YC

YC[1] <- X[1]*0.5 + rnorm(1, mean=0, sd=1)
YC

YC[2] <- X[2]*0.5 + (YC[1]-(X[1]*0.5)) + rnorm(1, mean=0, sd=1)
YC

for(i in 3:length(s))
	{
	YC[i] <- X[i]*0.5 + (YC[i-1]-(X[i-1]*0.5)) + rnorm(1, mean=0, sd=1)
	}
YC
plot(s, YC)

#1000 realizations of the process
for (k in 1:1000)
	{
	YC <- rep(NA, times=length(s))
	YC

	YC[1] <- X[1]*0.5 + rnorm(1, mean=0, sd=1)
	#YC

	YC[2] <- X[2]*0.5 + (YC[1]-(X[1]*0.5)) + rnorm(1, mean=0, sd=1)
	#YC

	for(i in 3:length(s))
		{
		YC[i] <- X[i]*0.5 + (YC[i-1]-(X[i-1]*0.5)) + rnorm(1, mean=0, sd=1)
		}

	#YC
	plot(s, YC, main = paste("Realization", k, " of the process"))
	Sys.sleep(1)
	}

#capture data generated by 1000 realizations of the process at selected spatial points 
yC.20 <- rep(NA, 1000)
yC.21 <- rep(NA, 1000)
yC.40 <- rep(NA, 1000)
yC.41 <- rep(NA, 1000)
yC.60 <- rep(NA, 1000)
yC.61 <- rep(NA, 1000)
yC.80 <- rep(NA, 1000)
yC.81 <- rep(NA, 1000)
 
for (k in 1:1000)
	{
	YC <- rep(NA, times=length(s))
	YC

	YC[1] <- X[1]*0.5 + rnorm(1, mean=0, sd=1)
	#YC

	YC[2] <- X[2]*0.5 + (YC[1]-(X[1]*0.5)) + rnorm(1, mean=0, sd=1)
	#YC

	for(i in 3:length(s))
		{
		YC[i] <- X[i]*0.5 + (YC[i-1]-(X[i-1]*0.5)) + rnorm(1, mean=0, sd=1)
		}

	#YC
	plot(s, YC, main = paste("Realization", k, " of the process"))
	#Sys.sleep(1)

	yC.20[k] <- YC[20]
	yC.21[k] <- YC[21]
	yC.40[k] <- YC[40]
	yC.41[k] <- YC[41]
	yC.60[k] <- YC[60]
	yC.61[k] <- YC[61]
	yC.80[k] <- YC[80]
	yC.81[k] <- YC[81]
	}

#examine the mean of the process at selected spatial points (s.20, s.21, ..., s.81) 
selected.s <- list(yC.20, yC.21, yC.40, yC.41, yC.60, yC.61, yC.80, yC.81)
boxplot(selected.s, names=c("s.20", "s.21", "s.40", "s.41", "s.60", "s.61", "s.80", "s.81"), ylab="Y")

#examine the covariance of the process between selected pairs of spatial points
plot(yC.20, yC.21)
cov(yC.20, yC.21)
cor.test(yC.20, yC.21)

plot(yC.40, yC.41)
cov(yC.40, yC.41)
cor.test(yC.40, yC.41)

plot(yC.60, yC.61)
cov(yC.60, yC.61)
cor.test(yC.60, yC.61)

plot(yC.80, yC.81)
cov(yC.80, yC.81)
cor.test(yC.80, yC.81)

plot(yC.20, yC.40)
cov(yC.20, yC.40)
cor.test(yC.20, yC.40)

plot(yC.20, yC.60)
cov(yC.20, yC.60)
cor.test(yC.20, yC.60)

plot(yC.20, yC.80)
cov(yC.20, yC.80)
cor.test(yC.20, yC.80)

plot(yC.40, yC.60)
cov(yC.40, yC.60)
cor.test(yC.40, yC.60)

plot(yC.40, yC.80)
cov(yC.40, yC.80)
cor.test(yC.40, yC.80)

plot(yC.60, yC.80)
cov(yC.60, yC.80)
cor.test(yC.60, yC.80)