# Load package
library(vars)


# Generate data
A = cbind(c(0.3,0.3),c(0,0.4))
B = cbind(c(-0.1,0.8),c(0,-0.3))

varstep <- function(A,B,x,y) {
  e = rnorm(2,0,1)
  A%*%x +B%*%y + e
}

x1 = c(0,0) 
y2 = c(0,0)

set.seed(1234)

results = cbind(y2,x1)
for (t in seq(1,50))
{
  temp <- x1
  x1 <- varstep(A,B,x1,y2)
  results <- cbind(results, x1)
  y2 <- temp
}

advertising = (results[1,1:50]*15+100)
sales = (results[2,1:50]*10+200)
plot(1:50, advertising, type = "l", ylim=c(0,250), xlab="Day", ylab="Advertising/Sales")
lines(1:50, sales, col="red")
legend(0, 65, legend=c("Sales", "Advertising"), col=c("red", "black"), lty=1, cex=0.8)

data<-cbind(advertising,sales)

# Clean up
rm(advertising)
rm(A)
rm(B)
rm(results)
rm(sales)
rm(t)
rm(temp)
rm(varstep)
rm(x1)
rm(y2)


# Estimate VAR model with 2 lags
model <- VAR(data, p = 2, type = "const")
summary(model)


# Get Impulse Response Functions from the model
feir1 <- irf(model, n.ahead = 10, ortho = FALSE, runs = 1000)

plot(feir1)
feir1

