.Rhistory

##graph the expected claims development pattern
plot(100*(rev(1/cumprod(rev(c(f, tail[tail>1.0001]))))), t="b",
main="Expected claims development pattern",
xlab="Dev. period", ylab="Development % of ultimate loss")
##fill in the remainder of the triangle
RAA <- cauto.tri
f <- c(f, f.tail)
fullRAA <- cbind(RAA, Ult = rep(0, 10))
for(k in 1:n){
fullRAA[(n-k+1):n, k+1] <- fullRAA[(n-k+1):n,k]*f[k]
}
round(fullRAA)
sum(fullRAA[ ,11] - getLatestCumulative(RAA))
linkratios <- c(attr(ata(RAA), "vwtd"), tail = 1.05)
round(linkratios, 3) # display to only three decimal places
LDF <- rev(cumprod(rev(linkratios)))
names(LDF) <- colnames(RAA) # so the display matches the triangle
round(LDF, 3)
currentEval <- getLatestCumulative(RAA)
# Reverse the LDFs so the first, least mature factor [1]
#	is applied to the last origin year (1990)
EstdUlt <- currentEval * rev(LDF) #
# Start with the body of the exhibit
Exhibit <- data.frame(currentEval, LDF = round(rev(LDF), 3), EstdUlt)
# Tack on a Total row
Exhibit <- rbind(Exhibit,
data.frame(currentEval=sum(currentEval), LDF=NA, EstdUlt=sum(EstdUlt),
row.names = "Total"))
Exhibit
lmCL <- function(i, Triangle){
lm(y~x+0, weights=1/Triangle[,i],
data=data.frame(x=Triangle[,i], y=Triangle[,i+1]))
}
sapply(lapply(c(1:(n-1)), lmCL, RAA), coef)
mack <- MackChainLadder(RAA, est.sigma="Mack")
mack
mack$f
mack$FullTriangle
plot(mack)
plot(mack, lattice=TRUE)
##Bootstrap
## See also the example in section 8 of England & Verrall (2002)
## on page 55.
B <- BootChainLadder(RAA, R=999, process.distr="gamma")
B
plot(B)
quantile(B, c(0.75,0.95,0.99, 0.995))
## fit a distribution to the IBNR
library(MASS)
plot(ecdf(B$IBNR.Totals))
## fit a log-normal distribution
fit <- fitdistr(B$IBNR.Totals[B$IBNR.Totals>0], "lognormal")
fit
curve(plnorm(x,fit$estimate["meanlog"], fit$estimate["sdlog"]),
col="red", add=TRUE)
##Clark Loss Development
ClarkLDF(RAA)
#calculate loss development after 20 years, or 240 months
ClarkLDF(RAA, maxage = 240)
ClarkLDF(RAA, G="weibull")
plot(ClarkLDF(RAA, G="weibull"))
ClarkCapeCod(RAA, Premium = 40000, G = "weibull")
plot(ClarkCapeCod(RAA, Premium = 40000, G = "weibull"))
# load data
#data(GenIns)
#GenIns <- GenIns / 1000
# fit Poisson GLM
RAA2 <- RAA/1000
(fit1 <- glmReserve(data=RAA))
##read in SNL commercial auto paid loss data for peer companies
cauto <- read.csv(file="/Users/robertmcpherson/Documents/workspace/ComAutoPeerCos_PaidLossSNL.csv")
options(prompt = "R> ", digits = 4, show.signif.stars = TRUE)
options(continue="   ")
lattice.options(default.theme = standard.theme(color = FALSE))
cauto.tri <- as.triangle(cauto,origin="AY",dev="Mos","RL")
cauto.tri
n <- 10
f <- sapply(1:(n-1),
function(i){
sum(cauto.tri[c(1:(n-i)),i+1])/sum(cauto.tri[c(1:(n-i)),i])
}
)
f
##graph the age-to-age factors
dev.period <- 1:(n-1)
plot(log(f-1) ~ dev.period, main="Log-linear extrapolation of age-to-age factors")
tail.model <- lm(log(f-1) ~ dev.period)
abline(tail.model)
co <- coef(tail.model)
## extrapolate another 100 dev. period
tail <- exp(co[1] + c((n + 1):(n + 100)) * co[2]) + 1
f.tail <- prod(tail)
f.tail
##graph the expected claims development pattern
plot(100*(rev(1/cumprod(rev(c(f, tail[tail>1.0001]))))), t="b",
main="Expected claims development pattern",
xlab="Dev. period", ylab="Development % of ultimate loss")
##fill in the remainder of the triangle
RAA <- cauto.tri
f <- c(f, f.tail)
fullRAA <- cbind(RAA, Ult = rep(0, 10))
for(k in 1:n){
fullRAA[(n-k+1):n, k+1] <- fullRAA[(n-k+1):n,k]*f[k]
}
round(fullRAA)
sum(fullRAA[ ,11] - getLatestCumulative(RAA))
linkratios <- c(attr(ata(RAA), "vwtd"), tail = 1.05)
round(linkratios, 3) # display to only three decimal places
LDF <- rev(cumprod(rev(linkratios)))
names(LDF) <- colnames(RAA) # so the display matches the triangle
round(LDF, 3)
currentEval <- getLatestCumulative(RAA)
# Reverse the LDFs so the first, least mature factor [1]
#	is applied to the last origin year (1990)
EstdUlt <- currentEval * rev(LDF) #
# Start with the body of the exhibit
Exhibit <- data.frame(currentEval, LDF = round(rev(LDF), 3), EstdUlt)
# Tack on a Total row
Exhibit <- rbind(Exhibit,
data.frame(currentEval=sum(currentEval), LDF=NA, EstdUlt=sum(EstdUlt),
row.names = "Total"))
Exhibit
lmCL <- function(i, Triangle){
lm(y~x+0, weights=1/Triangle[,i],
data=data.frame(x=Triangle[,i], y=Triangle[,i+1]))
}
sapply(lapply(c(1:(n-1)), lmCL, RAA), coef)
mack <- MackChainLadder(RAA, est.sigma="Mack")
mack
mack$f
mack$FullTriangle
plot(mack)
plot(mack, lattice=TRUE)
## See also the example in section 8 of England & Verrall (2002)
## on page 55.
B <- BootChainLadder(RAA, R=999, process.distr="gamma")
B
plot(B)
quantile(B, c(0.75,0.95,0.99, 0.995))
## fit a distribution to the IBNR
library(MASS)
plot(ecdf(B$IBNR.Totals))
## fit a log-normal distribution
fit <- fitdistr(B$IBNR.Totals[B$IBNR.Totals>0], "lognormal")
fit
curve(plnorm(x,fit$estimate["meanlog"], fit$estimate["sdlog"]),
col="red", add=TRUE)
##Clark Loss Development
ClarkLDF(RAA)
#calculate loss development after 20 years, or 240 months
ClarkLDF(RAA, maxage = 240)
ClarkLDF(RAA, G="weibull")
plot(ClarkLDF(RAA, G="weibull"))
ClarkCapeCod(RAA, Premium = 40000, G = "weibull")
plot(ClarkCapeCod(RAA, Premium = 40000, G = "weibull"))
# load data
#data(GenIns)
#GenIns <- GenIns / 1000
# fit Poisson GLM
RAA2 <- RAA/1000
(fit1 <- glmReserve(data=RAA))
ClarkLDF(RAA, G="weibull")
plot(ClarkLDF(RAA, G="weibull"))
ClarkCapeCod(RAA, Premium = 40000, G = "weibull")
plot(ClarkCapeCod(RAA, Premium = 40000, G = "weibull"))
# load data
#data(GenIns)
#GenIns <- GenIns / 1000
# fit Poisson GLM
RAA2 <- RAA/1000
(fit1 <- glmReserve(data=RAA))
RAA
(fit1 <- glmReserve(triangle=RAA))
str(RAA)
str(GenIns)
summary(fit1, type = "model")
###################################################
### code chunk number 69: ChainLadder.Rnw:1240-1244
###################################################
# Gamma GLM
(fit2 <- glmReserve(GenIns, var.power = 2))
# compound Poisson GLM (variance function estimated from the data):
#(fit3 <- glmReserve(GenIns, var.power = NULL))
###################################################
### code chunk number 70: ChainLadder.Rnw:1252-1254
###################################################
set.seed(11)
(fit5 <- glmReserve(GenIns, mse.method = "boot"))
###################################################
### code chunk number 71: ChainLadder.Rnw:1261-1262
###################################################
names(fit5)
###################################################
### code chunk number 72: ChainLadder.Rnw:1267-1271
###################################################
pr <- as.data.frame(fit5$sims.reserve.pred)
qv <- c(0.025, 0.25, 0.5, 0.75, 0.975)
res.q <- t(apply(pr, 2, quantile, qv))
print(format(round(res.q), big.mark = ","), quote = FALSE)
install.packages("Quandl")
install.packages("carat")
install.packages("VGAM")
install.packages("twilio")
install.packages("tseries")
rm(list=ls())
#setwd("R:/AnalyticsTeam/Personal/May/BitBucket/Gemini Time Series")
#setwd("C:/Users/rmcpherson/Documents/Segments/Phil Welt Segment/Gemini")
library(sqldf) #for running sql on data frames
library(dummies) #for creating one-hot encoding
library(forecast) #for the Holt-Winters forecast filter
library(glmnet) #for running regularized GLM
library(knitr) #for reproducible research, i.e., Markdown
#library(testthat)
library(BigVAR)
library(orderedLasso)
library(reshape)
library(ggplot2)
library(Quandl)
library(HDeconometrics)
library(imputeTS)
library(quantmod)
library(xts)
## ----set_globals---------------------------------------------------------
##########################
##Input Global Variables##
##########################
##########################
##Set the working directory
setwd("/Users/robertmcpherson/Documents/Code/R/TimeSeriesAnalysis/TimeSeries")
##########################
##########################
#Input the column name of the dependent variable to predict.
dependent.variable <- "DJI"
##########################
##########################
#Set the maximum lag for adjusting the variables in the data.
#each variable will get a new column for each lag, up to the maximum set here.
maxlag <- 24
##########################
##########################
#Type 'TRUE' if you want to include an offset in the GLM calculation, FALSE otherwise.
include.offset <- FALSE
##########################
##########################
#Frequency intervals: e.g., Yearly, Quarterly, Monthly, Daily
freq <- "quarterly"
##########################
##########################
#Input the column name that has the time increments in it, such as years, or year/months.
time.increment.variable <- "Date"
##########################
##########################
#Select whether to include plots with the arima, pre-whitening step
include.arima.plots <- FALSE
##########################
##########################
#Select whether to include cross correlation plots
include.cross.correlation.plots <- TRUE
##########################
##########################
#Select whether to include quartile to quartile (QQ) plots
include.QQ.plots <- FALSE
##########################
## ----load_data, results='hide'-------------------------------------------
#Note: this process takes the data in descending order, with the most recent data at the
CPI <- Quandl("RATEINF/CPI_USA", api_key="DJGcfzQc5RYP1JSycMBv", collapse=freq, start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)
CPI.lgth <- length(CPI[,1])
#Nonfinancial corporate business; short-term debt as a percentage of total debt, Annual
#shortTermDebtToLongTerm <- Quandl("FED/FL104140006_A", api_key="DJGcfzQc5RYP1JSycMBv", collapse=freq, start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)
#Financial soundness indicator, households; debt as a percent of gross domestic product
debtToGDP <- Quandl("FED/FL010000336_Q", api_key="DJGcfzQc5RYP1JSycMBv",collapse=freq, start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)
debtToGDP.lgth <- length(debtToGDP[,1])
#GDP <- Quandl("FED/FU086902001_A", api_key="DJGcfzQc5RYP1JSycMBv", collapse=freq, start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)
m1Velocity <- Quandl("FRED/M1V", api_key="DJGcfzQc5RYP1JSycMBv",collapse=freq, start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)
m1Velocity.lgth <- length(m1Velocity[,1])
m2Velocity <- Quandl("FRED/M2V", api_key="DJGcfzQc5RYP1JSycMBv", collapse=freq, start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)
m2Velocity.lgth <- length(m2Velocity[,1])
naturalUnemp <- Quandl("FRED/NROUST", api_key="DJGcfzQc5RYP1JSycMBv", collapse=freq, start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)
naturalUnemp.lgth <- length(naturalUnemp[,1])
realPotentialGDP <- Quandl("FRED/GDPPOT", api_key="DJGcfzQc5RYP1JSycMBv", collapse=freq, start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)
realPotentialGDP.lgth <- length(realPotentialGDP[,1])
totalPopulation <- Quandl("FRED/POP", api_key="DJGcfzQc5RYP1JSycMBv", collapse=freq, start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)
totalPopulation.lgth <- length(totalPopulation[,1])
#equitiesUncertaintyIndex <- Quandl("FRED/WLEMUINDXD", api_key="DJGcfzQc5RYP1JSycMBv", collapse=freq, start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)
#tradeWeightedUSDIndex <- Quandl("FRED/TWEXO", api_key="DJGcfzQc5RYP1JSycMBv", collapse=freq, start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)
#fiveYrForwardExpectedInflRate <- Quandl("FRED/T5YIFR", api_key="DJGcfzQc5RYP1JSycMBv", collapse=freq, start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)
#treasInflIndexedLongTermAvgYld <- Quandl("FRED/DLTIIT", api_key="DJGcfzQc5RYP1JSycMBv", collapse=freq, start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)
oneYrTreasSecondaryMkt <- Quandl("FRED/DTB1YR", api_key="DJGcfzQc5RYP1JSycMBv", collapse=freq, start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)
oneYrTreasSecondaryMkt.lgth <- length(oneYrTreasSecondaryMkt[,1])
x.min <- min(
CPI.lgth,
debtToGDP.lgth,
m1Velocity.lgth,
m2Velocity.lgth,
naturalUnemp.lgth,
realPotentialGDP.lgth,
totalPopulation.lgth,
oneYrTreasSecondaryMkt.lgth
)
CPI.trunc <- CPI[(CPI.lgth-x.min+1):CPI.lgth,]
debtToGDP.trunc <- debtToGDP[(debtToGDP.lgth-x.min+1):debtToGDP.lgth,]
m1Velocity.trunc <- m1Velocity[(m1Velocity.lgth-x.min+1):m1Velocity.lgth,]
m2Velocity.trunc <- m2Velocity[(m2Velocity.lgth-x.min+1):m2Velocity.lgth,]
naturalUnemp.trunc <- naturalUnemp[(naturalUnemp.lgth-x.min+1):naturalUnemp.lgth,]
realPotentialGDP.trunc <- realPotentialGDP[(realPotentialGDP.lgth-x.min+1):realPotentialGDP.lgth,]
totalPopulation.trunc <- totalPopulation[(totalPopulation.lgth-x.min+1):totalPopulation.lgth,]
oneYrTreasSecondaryMkt.trunc <- oneYrTreasSecondaryMkt[(oneYrTreasSecondaryMkt.lgth-x.min+1):oneYrTreasSecondaryMkt.lgth,]
raw.ts <- cbind(
CPI.trunc,
debtToGDP.trunc,
m1Velocity.trunc,
m2Velocity.trunc,
naturalUnemp.trunc,
realPotentialGDP.trunc,
totalPopulation.trunc,
oneYrTreasSecondaryMkt.trunc
)
#save time increment vector
#time.increments <- unique(raw_data_dummies[,time.increment.variable])
#time.increments <- time.increments[sort.list(time.increments, decreasing=FALSE)]
time.increments <- raw.ts[,time.increment.variable]
#rownames(raw.ts) <- raw.ts[,1]
data <- raw.ts[, !(names(raw.ts) %in% "Date")]
colnames(data) <- c("CPI","debtToGDP","m1Velocity","m2Velocity","naturalUnemp","realPotentialGDP","totalPopulation","oneYrTreasSecondaryMkt")
##Import birth data
#births <- read.csv("BirthsModified.csv")
#impute the na values
#births.interp <- na.interpolation(ts(births))
#births.export <- apply(births.interp, 2, rev)
#head(births.export)
#write.csv(births.export, "births_interp_no_blanks.csv")
#births.trunc <- births.interp[1:57,3:5]
#births.2009 <- apply(births.trunc, 2, rev)
#head(births.2009)
#na.interpolation(ts(births.trunc[,2]))
#births.interp <- apply(births.trunc, 2, na.interpolation)
#colnames(births.interp) <- c("Year", "Births", "BirthRate")
#head(births.interp)
#head(data)
#str(data)
#?rev
#births.2009 <- apply(births.interp[], 2, rev)
#################################
##Not included - series too short
##Get SPX data
#SPX <- getSymbols("^GSPC",auto.assign = FALSE)
#?getSymbols
#SPX.xts <- as.xts(SPX)
#SPX.yearly <- to.yearly(SPX.xts)
#getSymbols("DJI", src = "yahoo", from = start_date, to = end_date)
#DJI <- getSymbols("DJI", src = "yahoo")
#DJI.xts <- as.xts(DJI)
#DJI.interp <- na.interpolation(DJI.xts)
#DJI.yearly <- to.yearly(DJI.interp)
#################################
##Import historical Dow Jones Industrial Average data
#Dow <- read.csv("DJI_Historical_Chg_Data.csv")
#head(Dow)
#DowInterp <- na.interpolation(Dow)
#DJIChg <- Dow[1:57,2]
# Import stock data from Yahoo
start = as.Date("1960-12-30")
end = as.Date("2017-12-30")
getSymbols("DJI", src="yahoo", start=start, end=end)
DJI.xts <- as.xts(DJI)
DJI.interp <- na.interpolation(DJI.xts)
#DJI.yearly <- to.yearly(DJI.interp)
DJI.quarterly <- to.quarterly(DJI.interp)
DJI.close <- data.frame(DJI.quarterly$DJI.interp.Close)
colnames(DJI.close) <- "DJI"
# Adjust lengths of vectors to shortest length
index.a <- length(data[,1])
index.b <- length(DJI.close[,1])
data2 <- data.frame(data[((index.a - index.b)+1):index.a,])
#length(data2[,1])
#length(DJI.close)
#SeriesData <- cbind(births.2009, data, DJIChg)
#SeriesData <- cbind(births.2009, data, DJI)
SeriesData <- cbind(data2, DJI.close)
SeriesData
col.names <- colnames(SeriesData)
x <- SeriesData[, !(names(SeriesData) %in% dependent.variable)]
head(x)
#scale the independent variables
x.scaled <- scale(x)
x.scaled
x
scale(x)
x.scaled
#Isolate dependent variable values, based on name given in global variable inputs above
y <- SeriesData[,dependent.variable]
y.unscaled <- y
y.unscaled
#scale the dependent variable
y.scaled <- scale(y)
y.scaled
x.colnames <- colnames(x)
x.colnames
num.cols <- length(x[1,])
num.cols
x.arima.residuals = NULL
for (i in 1:num.cols){
fit <- auto.arima(x.scaled[,i])
pdf(file=paste("plots/", dependent.variable, "_arima_",x.colnames[i],".pdf",sep=""))
if(include.arima.plots == TRUE){
par(mar=c(8,4,2,2))
plot(forecast(fit,h=maxlag), sub=paste(x.colnames[i]))
}
dev.off()
#assemble a table of A
x.arima.residuals = NULL
for (i in 1:num.cols){
fit <- auto.arima(x.scaled[,i])
pdf(file=paste("plots/", dependent.variable, "_arima_",x.colnames[i],".pdf",sep=""))
if(include.arima.plots == TRUE){
par(mar=c(8,4,2,2))
plot(forecast(fit,h=maxlag), sub=paste(x.colnames[i]))
}
dev.off()
#assemble a table of ARIMA residuals for use in cross-correlation analysis
temp.resid <- resid(fit)
x.arima.residuals = NULL
for (i in 1:num.cols){
fit <- auto.arima(x.scaled[,i])
pdf(file=paste("plots/", dependent.variable, "_arima_",x.colnames[i],".pdf",sep=""))
if(include.arima.plots == TRUE){
par(mar=c(8,4,2,2))
plot(forecast(fit,h=maxlag), sub=paste(x.colnames[i]))
}
dev.off()
#assemble a table of ARIMA residuals for use in cross-correlation analysis
temp.resid <- resid(fit)
x.arima.residuals <- as.matrix(cbind(x.arima.residuals, temp.resid))
}
}
}
x.arima.residuals = NULL
for (i in 1:num.cols){
fit <- auto.arima(x.scaled[,i])
pdf(file=paste("plots/", dependent.variable, "_arima_",x.colnames[i],".pdf",sep=""))
if(include.arima.plots == TRUE){
par(mar=c(8,4,2,2))
plot(forecast(fit,h=maxlag), sub=paste(x.colnames[i]))
}
dev.off()
#assemble a table of ARIMA residuals for use in cross-correlation analysis
temp.resid <- resid(fit)
x.arima.residuals <- as.matrix(cbind(x.arima.residuals, temp.resid))
}
#run arima transformation on the dependent variable
fit=NULL
fit <- auto.arima(y.scaled)
par(mar=c(8,4,2,2))
pdf(file=paste("plots/",dependent.variable,"_arima.pdf",sep=""))
plot(forecast(fit,h=5), sub=paste(dependent.variable, sep=""))
dev.off()
y.arima.residuals <- resid(fit)
y.arima.residuals
if(include.QQ.plots == TRUE){
#check distributions of independent variables for normality
for (i in 1:length(x.scaled[1,])){
pdf(file=paste("plots/", dependent.variable, "_qqnorm_",x.colnames[i],".pdf",sep=""))
qqnorm(x.arima.residuals[,i], main=paste(x.colnames[i]))
dev.off()
}
#check dependent variable for normality
pdf(file=paste("plots/qqnorm_",dependent.variable,".pdf",sep=""))
qqnorm(y.arima.residuals, main=paste(dependent.variable,sep=""))
dev.off()
}
dir.create("plots")
if(include.cross.correlation.plots == TRUE){
for (i in 1:length(x[1,])){
pdf(file=paste("plots/", dependent.variable, "_ccf_",x.colnames[i],".pdf",sep=""))
par(mar=c(5,7,4,2)) #set the margins so title does not get cut off
ccf(x.arima.residuals[,i], y.arima.residuals, plot=TRUE, main=paste(x.colnames[i]), na.action = na.contiguous)
dev.off()
}
}
Y <- cbind.data.frame(x, y)
colnames(Y) <- col.names
nms <- dependent.variable
Y <- as.matrix(Y)
Y
# Fit a Basic VAR-L(3,4) on simulated data
T1=floor(nrow(Y)/3)
T2=floor(2*nrow(Y)/3)
T2
T1
# Break data into in and out of sample to test model accuracy
Yin = Y[1:T2,]
Yout = Y[(T2+1):(T1+T2),]
# BVAR
#?lbvar
#?predict
modelbvar=lbvar(Yin, p = 4, delta = 0.5)
predbvar=predict(modelbvar,h=4)
# Forecasts of the volatility
k=paste(dependent.variable)
pdf(file=paste("plots/", dependent.variable, "_forecast.pdf",sep=""))
plot(c(Y[,k],predbvar[,k]),type="l", main=paste(dependent.variable))
lines(c(rep(NA,length(Y[,k])),predbvar[,k]))
abline(v=length(Y[,k]),lty=2,col=4)
#dev.off()
# = Overall percentual error = #
MAPEbvar=abs((Yout-predbvar)/Yout)*100
aux=apply(MAPEbvar,2,lines,col="lightskyblue1")
lines(rowMeans(MAPEbvar),lwd=3,col=4,type="b")
dev.off()
install.packages(c("BigVAR", "imputeTS", "Quandl", "reshape"))
install.packages(c("BigVAR", "imputeTS", "Quandl", "reshape"))
Sys.which("pdflatex")