#************************************************************ 
# Simulation  for VAR Model
# install.packages("vars")
#------------------------------------------------------------

library(vars) # Load package

# Reset random number generator for reasons of reproducability
set.seed(2345) 
# Generate sample
t <- 200 # Number of time series observations
k <- 2 # Number of endogenous variables
p <- 2 # Number of lags
A.1 <- matrix(c(-.3,.6,-.4,.5),k) # Coefficient matrix of lag 1
A.2 <- matrix(c(-.1,-.2,.1,.05),k) # Coefficient matrix of lag 2
A <- cbind(A.1,A.2) # Companion form of the coefficient matrices

y <- matrix(0,k,t+2*p) # Raw series with zeros
for (i in (p+1):(t+2*p)){ # Generate series with e ~ N(0,0.5)
  y[,i] <- A.1%*%y[,i-1] + A.2%*%y[,i-2] + rnorm(k,0,.5)
}

y <- ts(t(y[,-(1:p)])) # Convert to time series format
names <- c("V1","V2") # Rename variables
plot.ts(y) # Plot the series

VARselect(y, lag.max=8)

var.1 <- VAR(y,2,type="none") # Estimate the model
summary(var.1)
ir.1 <- irf(var.1,impulse="series.1",response="series.2",n.ahead = 20,ortho = FALSE)
ir.1 <- irf(var.1,impulse="series.2",response="series.1",n.ahead = 20,ortho = FALSE)
plot(ir.1)
plot(ir.2)

ir.1 <- irf(var.1,impulse="Series.1",response="Series.2",n.ahead = 20,ortho = FALSE,cumulative=TRUE)
ir.2 <- irf(var.1,impulse="Series.2",response="Series.2",n.ahead = 20,ortho = FALSE)
plot(ir.1)
plot(ir.2)

par(mfrow=c(1, 2))
plot(ir.1)
plot(ir.2)

var.2 <- VAR(series,lag.max=5,type="none", ic="AIC") # Estimate the model
summary(var.2)
round(rbind(coef(var.1)[[1]][,1],coef(var.1)[[2]][,1]),2) # Rounded estimates