#################

G11 <- 1
G22 <- .5
sig <- sqrt(5)
rho <- -.25
G12 <- sqrt(G11*G22)*rho

correl <- matrix (NA, 10, 10)
diffs <- matrix(NA, 10, 10)
VAR <- rep(NA, 10)
    
for (j in 1:10){
	for (k in 1:10){
	    VAR[j] <- G11 + 2*G12*j + G22*j*j + sig^2
		correl[j,k] <- (G11 + G12*(j + k) + G22*j*k) / (sqrt(G11 + 2*G12*j + G22*j*j + sig^2)*
		                                           sqrt(G11 + 2*G12*k + G22*k*k + sig^2))
		diffs[j,k] <- abs(j-k)
	}
}
diffs

plot(c(1:10), VAR, main=paste("Variance"))
round(correl, 2)
par(mfrow=c(2,2))
plot(diffs[1,-1], correl[1,-1], main="t=1", ylim=c(0,1), xlab="DeltaT", ylab="rho(DeltaT)")
plot(diffs[2,-2], correl[2,-2], main="t=2", ylim=c(0,1), xlab="DeltaT", ylab="rho(DeltaT)")
plot(diffs[5,-5], correl[5,-5], main="t=5", ylim=c(0,1), xlab="DeltaT", ylab="rho(DeltaT)")
plot(diffs[9,-9], correl[9,-9], main="t=9", ylim=c(0,1), xlab="DeltaT", ylab="rho(DeltaT)")

par(mfrow=c(2,2))
plot(c(2:10), correl[1,-1], main="Cor(Y_{i1}, Y_{ij})", ylim=c(0,1), xlab="t_{ij}", ylab="rho(DeltaT)")
plot(c(1,3:10), correl[2,-2], main="Cor(Y_{i2}, Y_{ij})", ylim=c(0,1), xlab="t_{ij}", ylab="rho(DeltaT)")
plot(c(1:4,6:10), correl[5,-5], main="Cor(Y_{i5}, Y_{ij})", ylim=c(0,1), xlab="t_{ij}", ylab="rho(DeltaT)")
plot(c(1:8,10), correl[9,-9], main="Cor(Y_{i9}, Y_{ij})", ylim=c(0,1), xlab="t_{ij}", ylab="rho(DeltaT)")