## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
## ---- eval=FALSE--------------------------------------------------------------
# install_github("wbnicholson/BigVAR/BigVAR")
## -----------------------------------------------------------------------------
library(BigVAR)
data(Y)
## -----------------------------------------------------------------------------
# 3 x 7 coefficient matrix
B = BigVAR.fit(Y,struct='Basic',p=2,lambda=1)[,,1]
# construct 7 x 99 lag matrix of Y
Z = VARXLagCons(Y,p=2,oos=TRUE)$Z
# obtain out of sample forecasts
yhat = B%*%Z[,ncol(Z),drop=F]
## ---- echo=FALSE--------------------------------------------------------------
suppressMessages(library(BigVAR,quietly=TRUE,warn.conflicts=FALSE))
#onorm=c()
library(lattice)
#Oracle=diag(3)
#SigmaU=diag(3)
SparsityPlot <- function (B, p, k,s,m, title = NULL)
{
text <- c()
for (i in 1:p) {
text1 <- as.expression(bquote(bold(Phi)^(.(i))))
text <- append(text, text1)
}
## text <- c()
if(s>0){
for (i in (p+1):(p+s+1)) {
text1 <- as.expression(bquote(bold(beta)^(.(i-p))))
text <- append(text, text1)
}
}
f <- function(m) t(m)[, nrow(m):1]
rgb.palette <- colorRampPalette(c("white", "grey" ),space = "Lab")
## rgb.palette <- colorRampPalette(c("white", "blue"), space = "Lab")
at <- seq(k/2 + 0.5, p * (k)+ 0.5, by = k)
at2 <- seq(p*k+s/2+.5,p*k+s*m+.5,by=s)
at <- c(at,at2)
se2 = seq(1.75, by = k, length = k)
L2 <- levelplot(f(abs(B)), col.regions = rgb.palette, colorkey = NULL,
xlab = NULL, ylab = NULL, main = list(label = title,
cex = 1), panel = function(...) {
panel.levelplot(...)
panel.abline(a = NULL, b = 1, h = seq(1.5, m*s+p* k +
0.5, by = 1), v = seq(1.5, by = 1, length = p *
k+m*s))
bl1 <- seq(k + 0.5, p *
k + 0.5, by = k)
bl2 <- seq(p*k + 0.5, p *
k + 0.5+s*m, by = m)
b1 <- c(bl1,bl2)
panel.abline(a = NULL, b = 1, v = p*k+.5, lwd = 7)
panel.abline(a = NULL, b = 1, v = b1, lwd = 3)
}, scales = list(x = list(alternating = 1, labels = text,
cex = 1, at = at, tck = c(0, 0)), y = list(alternating = 0,
tck = c(0, 0))))
return(L2)
}
B1 <- matrix(rep(1,57)*rbinom(57,1,.6),nrow=3,ncol=19)
B2 <-matrix(0,nrow=3,ncol=19)
B2[,1:3] <- 1
B2[,10:12] <- 1
B2[,16] <- 1
B2[,19] <- 1
B3 <-matrix(0,nrow=3,ncol=19)
diag(B3[,1:3])<- 1
B3[,10:12] <- 1
diag(B3[,10:12])<- 0
B2[,16] <- 1
B2[,19] <- 1
B4 <-matrix(0,nrow=3,ncol=19)
B4[,1:3] <- rep(1,9)*rbinom(9,1,.6)
B4[,10:12] <- rep(1,9)*rbinom(9,1,.4)
B4[,18] <- c(1,0,1)
p=5;k=3
Lasso <- SparsityPlot(B1,p,k,m=2,s=2,title="Basic VARX-L")
Group <- SparsityPlot(B2,p,k,m=2,s=2,title="Lag VARX-L")
OO <- SparsityPlot(B3,p,k,m=2,s=2,title="Own/Other VARX-L")
Sparse <- SparsityPlot(B4,p,k,m=2,s=2,title="Sparse Lag VARX-L")
## ----echo=FALSE---------------------------------------------------------------
library(gridExtra,quietly=TRUE)
grid.arrange(Lasso,Group,OO,Sparse,ncol=2)
## ----echo=FALSE---------------------------------------------------------------
k=3;p=5
HLAGC <- matrix(0,nrow=k,ncol=k*p)
HLAGC[1,] <- 1
HLAGC[2,1:6] <- 1
HLAGC[3,1:12] <- 1
HLAGOO <- matrix(0,nrow=k,ncol=k*p)
HLAGOO[1,1:13] <- 1
HLAGOO[2,1:6] <- 1
HLAGOO[3,c(1:9,12)] <- 1
HLAGELEM <- matrix(0,nrow=k,ncol=k*p)
HLAGELEM[1,c(1:10,12,13)] <- 1
HLAGELEM[2,c(1,3,4,6,7,9,10,12,13,15)] <- 1
HLAGELEM[3,c(1:8,10:11,13:14)] <- 1
SparsityPlot <-
function (B, p, k,s,m, title = NULL)
{
text <- c()
for (i in 1:p) {
text1 <- as.expression(bquote(bold('\U03A6')^(.(i))))
text <- append(text, text1)
}
## text <- c()
if(m>0){
for (i in (p+1):(p+s+1)) {
text1 <- as.expression(bquote(bold('\U03B2')^(.(i-p))))
text <- append(text, text1)
}
}
f <- function(m) t(m)[, nrow(m):1]
rgb.palette <- colorRampPalette(c("white", "grey" ),space = "Lab")
## rgb.palette <- colorRampPalette(c("white", "blue"), space = "Lab")
at <- seq(k/2 + 0.5, p * (k)+ 0.5, by = k)
if(m>0){
at2 <- seq(p*k+m/2+.5,p*k+s*m+.5,by=m)}else{at2=c()}
at <- c(at,at2)
se2 = seq(1.75, by = k, length = k)
L2 <- levelplot(f(abs(B)), col.regions = rgb.palette, colorkey = NULL,
xlab = NULL, ylab = NULL, main = list(label = title,
cex = 1), panel = function(...) {
panel.levelplot(...)
panel.abline(a = NULL, b = 1, h = seq(1.5, m*s+p* k +
0.5, by = 1), v = seq(1.5, by = 1, length = p *
k+m*s))
bl1 <- seq(k + 0.5, p *
k + 0.5, by = k)
if(m>0){
bl2 <- seq(p*k + 0.5, p *
k + 0.5+s*m, by = m)}else(bl2=c())
b1 <- c(bl1,bl2)
panel.abline(a = NULL, b = 1, v = p*k+.5, lwd = 3)
panel.abline(a = NULL, b = 1, v = b1, lwd = 3)
}, scales = list(x = list(alternating = 1, labels = text,
cex = 1, at = at, tck = c(0, 0)), y = list(alternating = 0,
tck = c(0, 0))))
return(L2)
}
set.seed(1986)
B5 <-matrix(0,nrow=3,ncol=15)
B5[,1:3] <- rep(1,9)*rbinom(9,1,.85)
B5[,4:6] <- rep(1,9)*rbinom(9,1,.65)
B5[,7:9] <- rep(1,9)*rbinom(9,1,.45)
B5[,10:12] <- rep(1,9)*rbinom(9,1,.25)
B5[,13:15] <- rep(1,9)*rbinom(9,1,.05)
HV4 <- SparsityPlot(B5,5,k,0,0,title='Lag-Weighted Lasso')
HVC=SparsityPlot(HLAGC,p,k,0,0,title="Componentwise HLAG")
HVOO=SparsityPlot(HLAGOO,p,k,0,0,title="Own/Other HLAG")
HVELEM=SparsityPlot(HLAGELEM,p,k,0,0,title="Elementwise HLAG")
grid.arrange(HVC,HVOO,HVELEM,HV4,ncol=2)
## -----------------------------------------------------------------------------
data(Y)
# Create a Basic VAR-L (Lasso Penalty) with maximum lag order p=4, 10 grid points with lambda optimized according to rolling validation of 1-step ahead MSFE
mod1<-constructModel(Y,p=4,"Basic",gran=c(150,10),h=1,cv="Rolling",verbose=FALSE,IC=TRUE,model.controls=list(intercept=TRUE))
## -----------------------------------------------------------------------------
results=cv.BigVAR(mod1)
results
## -----------------------------------------------------------------------------
str(results)
## -----------------------------------------------------------------------------
plot(results)
## -----------------------------------------------------------------------------
mod2<-constructModel(Y,p=4,"Basic",gran=c(5,10),h=1,cv="Rolling",verbose=FALSE,IC=FALSE)
res2=cv.BigVAR(mod2)
plot(res2)
## -----------------------------------------------------------------------------
mod3<-constructModel(Y,p=4,"Basic",gran=c(500,10),h=1,cv="Rolling",verbose=FALSE,IC=FALSE)
res3=cv.BigVAR(mod3)
plot(res3)
## -----------------------------------------------------------------------------
SparsityPlot.BigVAR.results(results)
## -----------------------------------------------------------------------------
predict(results,n.ahead=1)
## -----------------------------------------------------------------------------
predict(results,n.ahead=1, confint=TRUE)
## -----------------------------------------------------------------------------
coef(results)
## ----echo=TRUE----------------------------------------------------------------
data(Y) # simulated multivariate time series
# coefficient matrix used to generate Y
data(Generator)
# note that coefficients with a darker shade are larger in magnitude
SparsityPlot(A[1:3,],p=4,3,s=0,m=0,title="Sparsity Structure of Generator Matrix")
## -----------------------------------------------------------------------------
# fit a Basic VARX-L with k=2,m=1,s=2,p=4,lambda=.01
VARX=list(k=2,s=2)
#returns k x (kp+ms+1) coefficient matrix
model=BigVAR.fit(Y,p,"Basic",lambda=1e-2,VARX=VARX,intercept=TRUE)
model
## -----------------------------------------------------------------------------
# N-fold cross validation for VAR
# Y: data
# nfolds: number of cross validation folds
# struct: penalty structure
# p: lag order
# nlambdas: number of lambdas:
# gran1: depth of lambda grid
# seed: set to make it reproducible
NFoldcv <- function(Y,nfolds,struct,p,nlambdas,gran1,seed)
{
MSFE <- matrix(0,nrow=nrow(Y),ncol=10)
A <- constructModel(Y,p,struct=struct,gran=c(gran1,nlambdas),verbose=F)
# construct lag matrix
Z1 <- VARXLagCons(Y,X=NULL,s=0,p=p,0,0)
trainZ <- Z1$Z[2:nrow(Z1$Z),]
trainY <- matrix(Y[(p+1):nrow(Y),],ncol=ncol(Y))
set.seed(seed)
inds <- sample(nrow(trainY))
B <- BigVAR.est(A)
lambda.grid <- B$lambda
folds <- cut(inds,breaks=nfolds,labels=FALSE)
MSFE <- matrix(0,nrow=nfolds,ncol=nlambdas)
for(i in 1:nfolds){
test <- trainY[which(folds==i),]
train <- trainY[which(folds!=i),]
testZ <-t(t(trainZ)[which(folds!=i),])
B=BigVAR.fit(train,p=p,lambda=lambda.grid,struct='Basic')
#iterate over lambdas
for(j in 1:nlambdas){
MSFETemp <- c()
for(k in 1:nrow(test)) {
tempZ <- testZ[,k,drop=FALSE]
bhat <- matrix(B[,2:dim(B)[2],j],nrow=ncol(Y),ncol=(p*ncol(Y)))
preds <- B[,1,j]+bhat%*%tempZ
MSFETemp <- c(MSFETemp,sum(abs(test[k,]-preds))^2)
}
MSFE[i,j] <- mean(MSFETemp)
}
}
return(list(MSFE=MSFE,lambdas=lambda.grid))
}
# 10 fold cv
MSFEs<-NFoldcv(Y,nfolds=10,"Basic",p=5,nlambdas=10,gran1=50,seed=2000)
# choose smaller lambda in case of ties (prevents extremely sparse solutions)
opt=MSFEs$lambda[max(which(colMeans(MSFEs$MSFE)==min(colMeans(MSFEs$MSFE))))]
opt
## -----------------------------------------------------------------------------
data(Y)
p <- 4
T1 <- floor(nrow(Y))/3
T2 <- floor(2*nrow(Y))/3
#Matrix of zeros for X
X <- matrix(0,nrow=nrow(Y),ncol=ncol(Y))
BICMSFE <- VARXForecastEval(Y,X,p,0,T1,T2,"BIC",1)
## -----------------------------------------------------------------------------
mod <- VARXFit(Y,3,NULL,NULL)
pred <-PredictVARX(mod)
pred
## -----------------------------------------------------------------------------
library(MASS)
k=3;p=6
B=matrix(0,nrow=k,ncol=p*k)
A1<- matrix(c(.4,-.07,.08,-.06,-.7,.07,-.08,.07,-.4),ncol=3,nrow=3)
A2 <- matrix(c(-.6,0,0,0,-.4,0,0,0,.5),ncol=3,nrow=3)
B[,1:k]=A1
B[,(5*k+1):(6*k)]=A2
A <- VarptoVar1MC(B,p,k)
set.seed(2000)
Y <-MultVarSim(k,A,p,.005*diag(k),500)
SparsityPlot(B,p,k,0,0, title='Sparsity Plot of VAR Coefficient Matrix')
## ----cache=TRUE---------------------------------------------------------------
library(MCS)
# train on first 250 observations
YTrain=Y[1:250,]
Loss <- c()
T1=1*floor(nrow(YTrain)/3);T2=2*floor(nrow(YTrain)/3)
p=8
structures<-c("Basic","BasicEN","Lag","SparseLag","OwnOther","HLAGC","HLAGOO","HLAGELEM","MCP","SCAD")
for(i in structures){
# construct BigVAR object; we will perform a dual grid search for the sparse lag and sparse own/other models
if(i%in%c("SparseLag","SparseOO")){
alpha=seq(0,1,length=10)
}else{
alpha=0
}
A<- constructModel(YTrain,p=p,struct=i,gran=c(100,10),T1=T1,T2=T2,verbose=FALSE,model.controls=list(intercept=FALSE,alpha=alpha))
# perform rolling cv
res<- cv.BigVAR(A)
# save out of sample loss for each structure
Loss <- cbind(Loss,res@OOSMSFE)
}
# construct AIC and BIC benchmarks
BIC <- VARXForecastEval(YTrain,matrix(0,nrow=nrow(YTrain)),p,0,T2,nrow(YTrain),"BIC",1)$MSFE
AIC <- VARXForecastEval(YTrain,matrix(0,nrow=nrow(YTrain)),p,0,T2,nrow(YTrain),"AIC",1)$MSFE
Loss <- as.data.frame(Loss)
names(Loss) <- structures
Loss <- cbind(Loss,BIC,AIC)
names(Loss)[(ncol(Loss)-1):ncol(Loss)] <- c("BIC","AIC")
names(Loss) <- paste0(names(Loss),"L")
mcs.test <- MCSprocedure(as.matrix(Loss),verbose=FALSE)
mcs.test
## -----------------------------------------------------------------------------
suppressMessages(library(expm))
# Phi k x kp coefficient matrix
# sigma kxk residual covariance matrix
# n number of time steps to run IRF
# p lag order
# k number of series
# Y0: k dimensional vector reflecting initialization of the IRF
generateIRF <- function(Phi,Sigma,n,k,p,Y0)
{
if(p>1){
A <-VarptoVar1MC(Phi,p,k)
}else{
A <- Phi
}
J <- matrix(0,nrow=k,ncol=k*p)
diag(J) <- 1
P <- t(chol(Sigma))
IRF <- matrix(0,nrow=k,ncol=n+1)
for(i in 0:n)
{
phi1 <- J%*%(A%^%i)%*%t(J)
theta20 <- phi1%*%P
IRF[,i+1] <- (theta20%*%Y0)
}
return(IRF)
}
## ---- echo=TRUE,eval=TRUE,cache=TRUE------------------------------------------
require(quantmod)
require(zoo)
# get GDP, Federal Funds Rate, CPI from FRED
#Gross Domestic Product (Relative to 2000)
getSymbols('GDP',src='FRED',type='xts')
GDP<- aggregate(GDP,as.yearqtr,mean)
GDP <- GDP/mean(GDP["2000"])*100
# Transformation Code: First Difference of Logged Variables
GDP <- diff(log(GDP))
index(GDP) <- as.yearqtr(index(GDP))
# Federal Funds Rate
getSymbols('FEDFUNDS',src='FRED',type='xts')
FFR <- aggregate(FEDFUNDS,as.yearqtr,mean)
# Transformation Code: First Difference
FFR <- diff(FFR)
# CPI ALL URBAN CONSUMERS, relative to 1983
getSymbols('CPIAUCSL',src='FRED',type='xts')
CPI <- aggregate(CPIAUCSL,as.yearqtr,mean)
CPI <- CPI/mean(CPI['1983'])*100
# Transformation code: difference of logged variables
CPI <- diff(log(CPI))
# Seasonally Adjusted M1
getSymbols('M1SL',src='FRED',type='xts')
M1<- aggregate(M1SL,as.yearqtr,mean)
# Transformation code, difference of logged variables
M1 <- diff(log(M1))
# combine series
Y <- cbind(CPI,FFR,GDP,M1)
names(Y) <- c("CPI","FFR","GDP","M1")
Y <- na.omit(Y)
k=ncol(Y)
T <- nrow(Y)
# start/end of rolling validation
T1 <- which(index(Y)=="1985 Q1")
T2 <- which(index(Y)=="2005 Q1")
#Demean
Y <- Y - (c(rep(1, nrow(Y))))%*%t(c(apply(Y[1:T1,], 2, mean)))
#Standarize Variance
for (i in 1:k) {
Y[, i] <- Y[, i]/apply(Y[1:T1,], 2, sd)[i]
}
library(expm)
# Fit an Elementwise HLAG model
Model1=constructModel(as.matrix(Y),p=4,struct="HLAGELEM",gran=c(25,10),verbose=FALSE,VARX=list(),T1=T1,T2=T2)
Model1Results=cv.BigVAR(Model1)
# generate IRF for 10 quarters following a 1 percent increase in the federal funds rate
IRFS <- generateIRF(Phi=Model1Results@betaPred[,2:ncol(Model1Results@betaPred)],Sigma=cov(Model1Results@resids),n=10,k=ncol(Y),p=4,Y0=c(0,.01,0,0))
IRFS <- generateIRF(Model1Results@betaPred[,2:ncol(Model1Results@betaPred)],cov(Model1Results@resids),10,4,4,c(0,.01,0,0))
## ----functions, include=FALSE, echo=FALSE-------------------------------------
# A function for captioning and referencing images
fig <- local({
i <- 0
ref <- list()
list(
cap=function(refName, text) {
i <<- i + 1
ref[[refName]] <<- i
paste("Figure ", i, ": ", text, sep="")
},
ref=function(refName) {
ref[[refName]]
})
})
## ----irf_plot,fig.cap='\\label{fig:irf_plot} Generated Impulse Reponses', echo=FALSE----
# Plot IRFs
par(mfrow=c(2,2))
for(i in 1:4)
{
plot(IRFS[i,],type='l',main=names(Y)[i],ylab="",xlab="")
}