example.R

# Simulating a VAR (Hoevenaars, 2005)
library(MASS)

a <- matrix(c(1.35/12,1.46/12, 7.04/12, -3.43/12, 5.11/12,
       1.21/12, 1.55/12), 7, 1)

B <- matrix(c(
c(0.34,0.01,0.00,-0.16,0.35,0.43,-0.56),
c(-0.18, -0.03, -0.07, -0.56, 1.20, 5.16, -1.82),
c(0.67, 0.03, 0.04, 5.25, -3.54, -0.54, 10.44),
c(-0.01, 0.00, 0.00, 0.95, 0.03, 0.02, -0.14),
c(-0.01, 0.01, 0.00, -0.01, 1.03, 0.19, -0.43),
c(0.02, -0.01, 0.00, 0.03, -0.07, 0.68, 0.54),
c(-0.01, 0.00, 0.00, -0.04, 0.05, 0.06, 0.78)),
7, 7, byrow=TRUE)

Sigma_ee <- 
matrix(c(
    c(0.58,0.40,0.27,-0.29,-0.35,0.17,0.14),
    c(0.40,4.60,0.20,-0.24,-0.66,0.15,0.61),
    c(0.27,0.20,7.57,-0.95,-0.07,-0.03,-0.07),
    c(-0.29,-0.24,-0.95,0.08,0.13,-0.03,0.09),
    c(-0.35,-0.66,-0.07,0.13,0.22,-0.83,-0.28),
    c(0.17,0.15,-0.03,-0.03,-0.83,0.17,-0.11),
    c(0.14,0.61,-0.07,0.09,-0.28,-0.11,0.07)),
    7,7,byrow=TRUE)

zt <- a

path <- matrix(0,7,1000+1)
path[,1] <- zt
for(i in 1:1000) {
    n <- Sigma_ee %*% mvrnorm(1,rep(0,7),Sigma = diag(7))
    zt <- B %*% zt + n 
    path[,i+1] <- zt
}

means <- c()
for(j in 1:7) {
    plot(ts(path[j,]))
    abline(h=a[j])
    means[j] <- mean(path[j,])
}

print(as.numeric(a))

## [1]  0.1125000  0.1216667  0.5866667 -0.2858333  0.4258333  0.1008333
## [7]  0.1291667

print(means)

## [1] -0.440259553  0.114609879  0.429371015 -0.972993879 -1.650826750
## [6]  0.260499428  0.004490577