# for plots use a large graphics window, R-3.1.1
# example for VAR
#setwd("C:/LVs/QF_FinancialEconometrics/R_Test")  

library(vars)

## Required
source("mq_res_R.txt")    

## read data  OR  use data of your own choice
data(Canada)           #   help(Canada)
#Canadian labor market, OECD data base,  1980Q1-2004Q4
# e    ... log of employment
# prod ... labor productivity 
# rw   ... real wages
# U    ... unemployment rate 

## a look at the data
n       <- nrow(Canada)
k       <- ncol(Canada)
# 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)
plot(Canada)


## estimate VAR ##
#help(VAR)

mod <- VAR(Canada, p = 2, type = "const", exogen = NULL )
summary(mod)    # roots are the modulus of the inverse roots of the lag polynomial 
coef(mod)       # coefficients
plot(mod)       # plot fitted, resids, ACF/PACF resids     
                # move on by enter/ clicking into the plot
                # make window large

# residuals
yf <- ts(fitted(mod), start= start(Canada)+c(0,2), frequency=frequency(Canada) ) #fitted / +p
ym <- ts(Canada[(mod$p+1):nrow(Canada), ],start= start(Canada)+c(0,2), frequency=frequency(Canada) )   #data / +p
ye <- ym - yf                               # residuals  
# OR   ye = resid(mod)
colnames(yf) <- colnames(ye) <- colnames(ym)
plot(ye)
mq_res(ye,16,mod$p)                             # mv Ljung-Box Test for residuals
# residual co-var-matrix
covar_e <- t(ye)%*%ye/(n - k*mod$p)
cor(ye)                                         # contemporary correlations

## AIC
k <- ncol(Canada);   p <- mod$p               # k , p
aic <- (-2)*logLik(mod)/(n-p) + 2*p*(k^2)/(n-p) 

# forecasting
plot( predict(mod, n.ahead= 6) )                
#fanchart( predict(mod, n.ahead= 6, ci= 0.95) )     

Phi(mod, nstep=3)      # matrices of the MA(infty) representation (of correlated errors)
                       # in our notation Psi (!)

## irf gives the the effect of 
# innovations of x with lag=0,1,2,...  to x_t,  to y_t, ...
# innovations of y with lag=0,1,2,...  to x_t,  to y_t, ...
irf(mod, n.ahead=8)          # impulse response fct orthogonalized (default) impulses
plot( irf(mod, n.ahead=6) )  # for move on enter


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


#print(mod)
#Psi(mod)               # in our notation Psi* (!)
                        # matrices in front of the orthogonalized errors

## EXERCISE:
#  (1) Set up a VAR (with data of your own choice).
#  (2) Choose the order of the VAR ( VAR(xy, ic ="AIC", ...) )
#  (3) Test the residuals wrt white noise.
#  (4) Comment on the roots of the characteristic polynomial.
#  (5) Generate forecasts for at least 4 periods into the future.
#  (6) Inspect the MA(infinity) representations of the model.
#  (7) Comment on the impulse response function.
#  (8) Comment on the forecast error decomposition.