#
# VAR - Simulated data
#
install.packages("vars")
install.packages("MASS")
library(MASS)
library(vars)

phi0 = c(-82,-193,187,12)
phi1=
matrix(c(
0.99 , 0.087   ,0.323 , -0.026,
0.145 ,  0.98  ,  0.524 ,  -0.079, 
-0.081, -0.193 ,  0.475 ,  0.075, 
-0.136, 0.064 ,-0.195 , 0.969),4,4,byrow=TRUE) 

Sigma = matrix(c( 
0.456 , 0.051 ,-0.028 , 0.025,
0.051 , 0.219 ,-0.134, -0.098,
-0.028, -0.134,  0.122 , 0.077,
0.025 ,-0.098 , 0.077 , 0.592),4,4)

set.seed(12354)
mu = rep(0,4)
n = 500
error = mvrnorm(2*n,mu,Sigma)
y = matrix(0,2*n,4)
for (t in 2:(2*n))
y[t,] = phi0+phi1%*%y[t-1,]+error[t,]

y = y[(n+1):(2*n),]

par(mfrow=c(2,2))
for (i in 1:4)
ts.plot(y[,i],ylab="",main=paste("Y",i,sep=""))


pairs(y)




# Choice of optimal lag length for an unrestricted VAR for a maximal lag length of eight. 
VARselect(y, lag.max = 8, type = "both")

p1ct = VAR(y,p=1) 

p1ct

parameters = rbind(p1ct$varresult$y1$coef,
p1ct$varresult$y2$coef,
p1ct$varresult$y3$coef,
p1ct$varresult$y4$coef)

phi0hat = as.vector(parameters[,5])
phi1hat = as.matrix(parameters[,1:4])

cbind(phi1,phi0)

plot(p1ct, names = "y1")

ser11 = serial.test(p1ct, lags.pt = 16, type = "PT.asymptotic") 
ser11$serial 

norm1 =  normality.test(p1ct) 
norm1$jb.mul


arch1 =  arch.test(p1ct, lags.multi = 5) 
arch1$arch.mul 

plot(arch1, names = "y1") 

plot(stability(p1ct), nc = 2)