# R-3.1.1, for plots use a large window
# exapmle for VAR

## packages ##
#install.packages("vars")
library(vars)

## Attention
source("mq_R.txt")    # in working dir / mv Ljung-Box Test / Call: mq(xx,16) 

## data ##
data(Canada)       #  read data;  help(Canada)
## a look at the data
# n       <- length(Canada[,1])
# d.dat  <- rbind( Canada[1:3,], Canada[(n-2):n,] )
# d.time <- c( index(Canada)[1:3], index(Canada)[(n-2):n] )
# d.info  <- cbind( d.time, d.dat)
# d.info
plot(Canada)

## generated data
#AR(1)
n <- 200; x1 <- c(1:n)*0; x1 <- as.ts(x1);
aa <- 0.5; x1[1]=0.1; for (i in (2:n)) { x1[i]=aa*x1[i-1] + rnorm(1)}; x1 <- as.ts(x1); 
## WN
#x1 <- as.ts(rnorm(n)); 
x2 <- as.ts(rnorm(n)); x3 <- as.ts(rnorm(n)); x4 <- as.ts(rnorm(n)); 
xx <- cbind(x1,x2,x3,x4)                   # data.class(xx) is "mts"
# plot(xx)
mq(xx,16) 


## estimate VAR ##
#help(VAR)

m <- VAR(Canada, p = 2, type = "const", exogen = NULL )
summary(m)            # roots are the modulus of the inverse roots of the lag polynomial 
                            # estimation method not given
                            # we loose the first p observations in each series
coef(m)                              # coefficients
plot(m)                               # plot fitted, resids, ACF/PACF resids: move on by enter

# residuals
yf <- ts(fitted(m), start= start(Canada)+c(0,2), frequency=frequency(Canada) )    #fitted / +p
ym <- ts(Canada,  start= start(Canada)+c(0,2), frequency=frequency(Canada) )    #data /  +p
ye <- ym - yf                               # residuals  /   res = ts( resid(m), ....)
colnames(yf) <- colnames(ye) <- colnames(ym)
plot(ye)
mq(ye,16)                                   # mv Ljung-Box Test

## AIC
k <- 4;   p <- 2               # k , p
n.m <- n - p
aic <- (-2)*logLik(m)/n.m + 2*p*(k^2)/n.m 

# forecasting
plot( predict(m, n.ahead= 8) )

Phi(m, nstep=3)  # matrices of the MA(infty) representation

irf(m, n.ahead=8)  # impulse response fct
plot( irf(m, n.ahead=8) )  # for move on enter

fevd(m, n.ahead=8)   # forecast error decomposition wrt orthogonalized residuals
plot( fevd(m, n.ahead=8) )




print(m)
Psi(m)