varmat <- matrix(c(1, .7, 0, 0,
                  .7, 1, 0, 0,
                  0, 0, 1, .7,
                  0, 0, .7, 1), ncol=4) #, byrow = T)
periods <- 100
#data <- rmvnorm(1000, mean = rep(0,4), 
#                cov = varmat)
#var(data)

data <- array(NA, dim=c(300, 4, periods))
data[,,1] <- rmvnorm(300, mean = rep(0,4), 
                cov = varmat)

for (j in 2:periods) {
	data[,1,j] <- (0.75 * data[,1,j-1]) + rnorm(300, 0, sd = 1)
	data[,2,j] <- (0.70 * data[,2,j-1]) + 0.2 * data[,1,j] 
	               #+ 10 # ifelse(j<=10, 0, 10)
	               + rnorm(300, 0, sd = 1)
	data[,3,j] <- (0.75 * data[,3,j-1]) + rnorm(300, 0, sd = 1)
	data[,4,j] <- (0.70 * data[,4,j-1]) + 0.1 * data[,1,j] 
	               #-10 # ifelse(j<=10, 0, -10)
	               + rnorm(300, 0, sd = 1)
}
x <-apply(data, 3, cor)
par(mfrow=c(2,1))
hist(x[2,])
hist(x[12,])
par(mfrow=c(2,1))
plot(c(1, periods), c(min(data[,2,]), max(data[,2,])), type="n")
for (i in 1:300) lines(1:periods, data[i,2,])

plot(c(1, periods), c(min(data[,4,]), max(data[,4,])), type="n")
for (i in 1:300) lines(1:periods, data[i,4,])


par(mfrow=c(4,2))
plot(data[1,1,])
plot(data[10,1,])
plot(data[1,2,])
plot(data[10,2,])
plot(data[1,3,])
plot(data[10,3,])
plot(data[1,4,])
plot(data[10,4,])

matlines(1:20, t(data[,i,]))


apply(data, 3, cor)