vars package of R - AIC after restrict

1.7k views Asked by At

I fit a vector autoregression model on the Canada data in the vars package and then restrict based on a t-value of 1.64.

library(vars)
data("Canada")
var.can1 <- VAR(Canada, p = 2, type = "none")

summary(var.can1)

VAR Estimation Results:
========================= 
Endogenous variables: e, prod, rw, U 
Deterministic variables: none 
Sample size: 82 
Log Likelihood: -184.045 
Roots of the characteristic polynomial:
    1 0.9783 0.9113 0.9113 0.7474 0.1613 0.1613 0.1572
Call:
VAR(y = Canada, p = 2, type = "none")

# AIC BIC etc.
VARSelect(Canada, lag.max = 2, type = "none")$criteria

var.can2 <- restrict(var.can1, method = "ser", thresh = 1.64)

summary(var.can2)

VAR Estimation Results:
========================= 
Endogenous variables: e, prod, rw, U 
Deterministic variables: none 
Sample size: 82 
Log Likelihood: -191.376 
Roots of the characteristic polynomial:
    1 0.9742 0.9272 0.9272 0.7753 0.2105 0.2105 0.005071
Call:
VAR(y = Canada, p = 2, type = "none")

I then want to obtain the revised information criteria, but cannot see a way to do this. Does anyone know how?

EDIT 1

So I try to derive the AIC for the unrestricted model:

vars::VARselect(Canada, lag.max = 2, type = "none")$criteria
                  1            2
AIC(n) -5.600280680 -6.082112784
HQ(n)  -5.411741957 -5.705035337
SC(n)  -5.130676924 -5.142905272
FPE(n)  0.003697972  0.002289041

s <- summary(var.can1)
s$covres
                e          prod          rw             U
e     0.140560073  0.0056629572 -0.03893668 -0.0798565366
prod  0.005662957  0.4358209615  0.06689687 -0.0005118419
rw   -0.038936678  0.0668968657  0.60125872  0.0309232731
U    -0.079856537 -0.0005118419  0.03092327  0.0899478736

From New introduction to multiple time series analysis Luetkepohl, Helmut 2007, pg 147:

$$AIC(m) = ln(det(covres)) + \frac{2mk^2}{T}$$

m is lag order, k is number of series, T is sample size

But I get:

-6.451984 + 2*2*4^2/84 = -5.69

which does not equal -5.600280680

1

There are 1 answers

0
t-student On BEST ANSWER

Digging around in the code I find that the residual covariance matrix that is reported in the summary is not what is actually used to compute the AIC.

Pretty frustrating, some would say a bug.

library(vars)
data("Canada")
var.can1 <- VAR(Canada, p = 2, type = "none")

s <- summary(var.can1)

# Variance covariance matrix for the resid
s$covres

# e          prod          rw             U
# e     0.140560073  0.0056629572 -0.03893668 -0.0798565366
# prod  0.005662957  0.4358209615  0.06689687 -0.0005118419
# rw   -0.038936678  0.0668968657  0.60125872  0.0309232731
# U    -0.079856537 -0.0005118419  0.03092327  0.0899478736

vars::VARselect(Canada, lag.max = 2, type = "none")$criteria


# AIC is defined as:
# AICm = ln(det(sigma)) + (2pk^2)/N
# p = lag order, K = num series

p <- 2
K <- 4

N <- nrow(Canada) - p
(AIC1 <- log(det(s$covres))  + (2*2*4^2)/N)

# [1] -5.671496 is nothing like -6.082112784


# The residual covariance matrix
# that is reported in the summary is not what is actually used
# to compute the AIC. 

myresid <- residuals(var.can1)
(mysigma <- crossprod(myresid) *  (1/N))

# e          prod          rw             U
# e     0.12684690  0.0051104797 -0.03513798 -0.0720656604
# prod  0.00511048  0.3933018507  0.06037034 -0.0004619127
# rw   -0.03513798  0.0603703437  0.54259933  0.0279063671
# U    -0.07206566 -0.0004619127  0.02790637  0.0811724773

log(det(mysigma)) + (2* p * K^2)/N
# [1] -6.082113, which is very like -6.082112784