1library(vars)
2
3data <- read.csv('/home/wesm/code/statsmodels/scikits/statsmodels/datasets/macrodata/macrodata.csv')
4names <- colnames(data)
5data <- log(data[c('realgdp', 'realcons', 'realinv')])
6data <- sapply(data, diff)
7
8reorder.coefs <- function(coefs) {
9 n <- dim(coefs)[1]
10 # put constant first...
11 coefs[c(n, seq(1:(n-1))),]
12}
13
14extract.mat <- function(lst, i) {
15  sapply(lst, function(x) x[,i])
16}
17
18get.coefs <- function(est) {
19  reorder.coefs(extract.mat(coef(est), 1))
20}
21
22get.stderr <- function(est) {
23  reorder.coefs(extract.mat(coef(est), 2))
24}
25
26reorder.phi <- function(phis) {
27  # Puts things in more proper C order for comparison purposes in Python
28
29  k <- dim(phis)[1]
30  n <- dim(phis)[3]
31
32  arr <- array(dim=c(n, k, k))
33
34  for (i in 1:n)
35    arr[i,,] <- phis[,,i]
36
37  arr
38}
39
40causality.matrix <- function(est) {
41  names <- colnames(est$y)
42  K <- est$K
43
44  # p-values
45  result <- matrix(0, nrow=K, ncol=)
46  for (i in 1:K) {
47    ## # causes
48    ## result[i,1] <- causality(est, cause=names[i])$Granger$p.value
49
50    # caused by others
51    result[i,1] <- causality(est, cause=names[-i])$Granger$p.value
52  }
53
54  colnames(result) <- c("causedby")
55
56  result
57}
58
59get.results <- function(data, p=1) {
60  sel <- VARselect(data, p) # do at most p
61
62  est <- VAR(data, p=p)
63
64  K <- ncol(data)
65
66  nirfs <- 5
67  orth.irf <- irf(est, n.ahead=nirfs, boot=F)$irf
68  irf <- irf(est, n.ahead=nirfs, boot=F, orth=F)$irf
69
70  crit <- t(sel$criteria)
71  colnames(crit) <- c('aic', 'hqic', 'sic', 'fpe')
72
73  resid <- resid(est)
74  detomega <- det(crossprod(resid) / (est$obs - K * p - 1))
75
76  n.ahead <- 5
77
78  list(coefs=get.coefs(est),
79       stderr=get.stderr(est),
80       obs=est$obs,
81       totobs=est$totobs,
82       type=est$type,
83       crit=as.list(crit[p,]),
84       nirfs=nirfs,
85       orthirf=orth.irf,
86       irf=irf,
87       causality=causality.matrix(est),
88       detomega=detomega,
89       loglike=as.numeric(logLik(est)),
90       nahead=n.ahead,
91       phis=Phi(est, n.ahead))
92}
93
94k <- dim(data)[2]
95result <- get.results(data, p=2)
96
97est = VAR(data, p=2)
98