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