LeadingIndicatorTimeSeries_Example.Rnw

\documentclass[12pt]{article}         % the type of document and font size (default 10pt)
\usepackage[margin=1.0in]{geometry}   % sets all margins to 1in, can be changed
\usepackage{moreverb}                 % for verbatimtabinput -- LaTeX environment
\usepackage{url}                      % for \url{} command
\usepackage{amssymb}                  % for many mathematical symbols
\usepackage[pdftex]{lscape}           % for landscaped tables
\usepackage{longtable}                % for tables that break over multiple pages
\usepackage{graphicx}
\title{Consumer Price Index Forecast with Demographic and Economic Leading Indicators}  % to specify title
\author{Bob McPherson}          % to specify author(s)
\begin{document}                      % document begins here
\SweaveOpts{concordance=TRUE}

% If .nw file contains graphs: To specify that EPS/PDF graph files are to be
% saved to 'graphics' sub-folder
%     NOTE: 'graphics' sub-folder must exist prior to Sweave step
%\SweaveOpts{prefix.string=graphics/plot}

% If .nw file contains graphs: to modify (shrink/enlarge} size of graphics
% file inserted
%         NOTE: can be specified/modified before any graph chunk
\setkeys{Gin}{width=1.0\textwidth}

\maketitle              % makes the title
%\tableofcontents        % inserts TOC (section, sub-section, etc numbers and titles)
%\listoftables           % inserts LOT (numbers and captions)
%\listoffigures          % inserts LOF (numbers and captions)
%                        %     NOTE: graph chunk must be wrapped with \begin{figure},
%                        %  \end{figure}, and \caption{}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}


<<echo=FALSE, results=hide>>=

#install.packages("sqldf")
#install.packages("dummies")
#install.packages("forecast")
#install.packages("orderedLasso")
#install.packages("glmnet")
#install.packages("glmnetcr")
#install.packages("h2o")
#install.packages("addendum")
#install.packages("testthat")
#devtools::use_testthat

rm(list=ls())

setwd("<<<PUT PATH HERE>>>")
dir.create("plots")

library(tidyr)
library(dplyr)
library(sqldf) #for running sql on data frames
library(dummies) #for creating one-hot encoding
library(forecast) #for the Holt-Winters forecast filter
#library(orderedLasso)
library(glmnet) #for running regularized GLM
library(glmnetcr) #for running regularized GLM
#library(h2o)
library(knitr) #for reproducible research, i.e., Markdown
library(testthat)
library(xtable)
library(HDeconometrics)
#ls("package:HDeconometrics")

library(MTS) #https://www.rdocumentation.org/packages/MTS/versions/1.0
#also: https://www.rdocumentation.org/packages/MTS/versions/1.0/topics/BVAR
#ls("package:MTS")
#?GrangerTest
#?BVAR

library(VARsignR) #https://cran.r-project.org/web/packages/VARsignR/vignettes/VARsignR-vignette.html

#also check out BMR and see if I can find the relevant package: https://www.kthohr.com/bmr/BMR.pdf



##########################
##Input Global Variables##
##########################

##########################
#Input the column name of the dependent variable to predict.
dependent.variable <- "total_incurred"
#dependent.variable <- "WP"
##########################

##########################
#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 <- 36
##########################

##########################
#Input the column name that has the time increments in it, such as years, or year/months.
time.increment.variable <- "policy_month"
##########################

##########################
#Select whether to include plots with the arima, pre-whitening step
include.arima.plots <- TRUE
##########################

##########################
#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
##########################

#Note: this process takes the data in descending order, with the most recent data at the
#bottom, or end of the list/table.

#dataset <- data.frame(read.csv(file="C:/Users/rmcpherson/Documents/Segments/Phil Welt Segment/Gemini/GeminiTimeSeries20170628.csv"))

load("NIC_property_data.rda")
str(dat)

#remove the most recent year of data, as it does not represent a complete year's worth of data
#dataset <- dataset[dataset$EffYear!= 2017,]

head(dat)

write.csv(head(dat), "sampledata.csv")

length(dat[,1])

sink("summary.txt")
summary(dat)
sink()

##One-Hot Encoding: Creating Dummy Variables

#turn categorical variables into dummy variables, and keep numerical data unchanged
raw_data_dummies <- dummy.data.frame(data=dat, names=c(
"policy_status"
,"form"
,"Package"
,"fdescription"
,"grp3description"
,"grp2description"
,"wrbc_asl"
,"protection_class"
,"construction"
,"AOP_Deductible"
,"location_building"
,"wind"
,"uw_scale"
,"renewal"
,"ProgramType"
,"ProgramDesc"
,"NationalAgency"
,"state"
#,"terr_code" #this variable has category 'C' also appearing as lower case - results in duplication error
,"MP.flag"
,"row_num"
), sep=".", all=TRUE)

#run again, but leave out all data types except the dummy variables, by changing to all=FALSE
#turn categorical variables into dummy variables, and keep numerical data unchanged
raw_data_dummies2 <- dummy.data.frame(data=dat, names=c(
"policy_status"
,"form"
,"Package"
,"fdescription"
,"grp3description"
,"grp2description"
,"wrbc_asl"
,"protection_class"
,"construction"
,"AOP_Deductible"
,"location_building"
,"wind"
,"uw_scale"
,"renewal"
,"ProgramType"
,"ProgramDesc"
,"NationalAgency"
,"state"
#,"terr_code" #this variable has category 'C' also appearing as lower case - results in duplication error
,"MP.flag"
,"row_num"
), sep=".", all=FALSE)


#head(raw_data_dummies)


#sum the ones (i.e., exposure counts) for each category by time period
cats <- aggregate(raw_data_dummies2, by=list(dat$policy_yr_qtr), FUN=sum)
head(cats)

#isolate all the numeric variables (i.e., not based on categorical data)
raw_data_numeric <- data.frame(dat[,c(
"WH_Deductible"
,"LOI"
,"Year_Built"
,"miles_to_ocean"
,"baserate"
,"LOI_Factor"
,"uw_scale_factor"
,"age_factor"
,"construction_factor"
,"coverage_factor"
,"deductible_factor"
,"form_factor"
,"protection_class_factor"
,"regionfactor"
,"windload"
,"minrate"
,"minimum_premium"
,"total_factor"
,"Charged.Prem"
,"Manual.Prem"
,"Policy.Premium"
,"loss_paid"
,"lae_paid"
,"loss_reserve"
,"lae_reserve"
,"total_incurred"
,"clm_ct"
,"capped_loss_ratio"
,"policy_yr_qtr"
)])

#aggregate numeric data by averaging over time periods
nums <- aggregate(raw_data_numeric, by=list(raw_data_numeric$policy_yr_qtr), FUN=mean)
#?apply
#head(nums)

#aggregate numeric data by summing over time periods
nums.sum <- aggregate(raw_data_numeric, by=list(raw_data_numeric$policy_yr_qtr), FUN=sum)

incurred.loss.ratio <- nums.sum$total_incurred/nums.sum$Charged.Prem
#clms.freq <- nums.sum$clm_ct

SeriesData <- cbind(cats, nums)
#head(SeriesData)
write.csv(SeriesData, "SeriesData.csv")

#fix column names to have proper name syntax
tidy.colnames <- make.names(colnames(SeriesData), unique=TRUE)
colnames(SeriesData) <- tidy.colnames

#get list of variables, and paste into exported list below
write.csv(file="colnames.csv",x=tidy.colnames)

#Use the list below as a starting point for selecting predictor variables. Uncomment variables to select.

x <- SeriesData


#head(x)

#scale the independent variables
x.scaled <- scale(x)

#Isolate dependent variable values, based on name given in global variable inputs above
y <- SeriesData[,dependent.variable]
y.unscaled <- y

#scale the dependent variable
y.scaled <- scale(y)

#save column names
x.colnames <- data.frame(colnames(x))

##ARIMA Time Series Analysis

#i=1
num.cols <- length(x[1,])
#apply(x,1,function(x) sum(is.na(x)))
#str(x)
#?auto.arima
#generate ARIMA plots...intent is to get ARIMA parameters, rather than forecasts
x.arima.residuals = NULL
for (i in 1:num.cols){
  fit <- auto.arima(x.scaled[,i])
  if(include.arima.plots == TRUE){
     pdf(paste("plots/ARIMA_",x.colnames[i,],".pdf", sep=""))
     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)
pdf(paste("plots/ARIMA_",dependent.variable,".pdf", sep=""))
par(mar=c(8,4,2,2))
plot(forecast(fit,h=1), sub=paste(dependent.variable, sep=""))
dev.off()
y.arima.residuals <- resid(fit)

#create a standardized, scaled, and normalized version of the data
#?scale
#glm

if(include.QQ.plots == TRUE){
#check distributions of independent variables for normality
  for (i in 1:length(x.scaled[1,])){
    pdf(paste("plots/QQ_",x.colnames[i,],".pdf", sep=""))
    qqnorm(x.arima.residuals[,i], main=paste(x.colnames[i,]))
    dev.off()
  }
}

#check dependent variable for normality
#qqnorm(y.arima.residuals, main=paste(dependent.variable,sep=""))

#check offset variable for normality
#qqnorm(offset.arima.residuals, main=paste(offset.variable,sep=""))

##Cross Correlation Analysis

#i=1
##cross correlation analysis
#leading indicators in 'x' will have negative lag values for the most significant
#correlations in the chart.
#note: analysis is run on ARIMA residuals so as to pre-whiten the data

##function for generating cross correlation tables and plots
cross.correl <- function(indep.vars.prewhitened, dep.vars.prewhitened, plots.subdirectory.name){
#rm(cross.correl)
  x <- indep.vars.prewhitened
  y <- dep.vars.prewhitened
  subdir <- plots.subdirectory.name
  dir.create(paste("plots/",subdir,sep=""))

  pos.cor.tbl <- NULL
  neg.cor.tbl <- NULL
  tmp <- NULL
  for (i in 1:length(x[1,])){
    cross.correl <- ccf(x[,i], y, plot=FALSE, na.action = na.contiguous)

    #find best correlation
    ind.max <- which(abs(cross.correl$acf[1:length(cross.correl$acf)])==max(abs(cross.correl$acf[1:length(cross.correl$acf)])))
    #extract optimal lag, and optimal corresponding correlation coefficient
    max.cor <- cross.correl$acf[ind.max]
    lag.opt <- cross.correl$lag[ind.max]

    #calculate statistical significance of the optimal correlation
    p.val <- 2 * (1 - pnorm(abs(max.cor), mean = 0, sd = 1/sqrt(cross.correl$n.used)))

    ## positively correlated, statistically significant, leading indicators
    if(p.val <= 0.05 && lag.opt < 0 && max.cor > 0){
       #make table
       tmp <- cbind(paste(x.colnames[i,]), round(max.cor,2), lag.opt, round(p.val,3))
       pos.cor.tbl <- rbind(tmp, pos.cor.tbl)
       #make plot
       pdf(paste("plots/",subdir,"/CCF_pos_",x.colnames[i,],".pdf", sep=""))
       par(mar=c(5,7,4,2)) #set the margins so title does not get cut off
       ccf(x[,i], y, plot=TRUE, main=paste(x.colnames[i,]), na.action = na.contiguous)
       dev.off()
    } #end if

    ## negatively correlated, statistically significant, leading indicators
    if(p.val <= 0.05 && lag.opt < 0 && max.cor < 0){
       #make table
       tmp <- cbind(paste(x.colnames[i,]), round(max.cor,2), lag.opt, round(p.val,3))
       neg.cor.tbl <- rbind(tmp, neg.cor.tbl)
       #make plot
       pdf(paste("plots/",subdir,"/CCF_neg_",x.colnames[i,],".pdf", sep=""))
       par(mar=c(5,7,4,2)) #set the margins so title does not get cut off
       ccf(x[,i], y, plot=TRUE, main=paste(x.colnames[i,]), na.action = na.contiguous)
       dev.off()
    } #end if
} #end loop

  ##export csv reports:
  #one for significant positive leading indicators, and one for significant negative leading indicators
  #positive correlation leading indicator summary
  colnames(pos.cor.tbl) <- c("Variable", "Cor", "Lag", "p_val")
  print(xtable(data.frame(pos.cor.tbl),caption="Positively correlated leading indicators", label="Table:positiveleading"))
  write.csv(data.frame(pos.cor.tbl), paste("plots/",subdir,"/LeadingIndicators_Positive.csv",sep=""))

  #negative correlation leading indicator summary
  colnames(neg.cor.tbl) <- c("Variable", "Cor", "Lag", "p-val")
  print(xtable(data.frame(neg.cor.tbl),caption="Negatively correlated leading indicators", label="Table:negativeleading"))
  write.csv(data.frame(neg.cor.tbl), paste("plots/",subdir,"/LeadingIndicators_Negative.csv",sep=""))

  #combine positive and negative leading indicator lists into one reference table
  leading.indicators <- rbind(pos.cor.tbl, neg.cor.tbl)
  return(leading.indicators)
} #end function

#make list of leading indicators with TOTAL CLAIMS as target variable - for reference - using above function
  leading.indicators.total.clms <- cross.correl(indep.vars.prewhitened = x.arima.residuals, dep.vars.prewhitened = y.arima.residuals, plots.subdirectory.name = "total_claims")

#make list of leading indicators with INCURRED LOSS RATIO as target variable - using above function
  leading.indicators.incurred.lr <- cross.correl(indep.vars.prewhitened = x.arima.residuals, dep.vars.prewhitened = incurred.loss.ratio, plots.subdirectory.name = "incurred_loss_ratio")

#write list of leading indicators to csv file
write.csv(leading.indicators.incurred.lr, "list.leading.indicators.incurred.lr.csv")

##assemble matrix of all significant leading indicators and corresponding data
lead.ind.inc.lr.matrix <- x.scaled[,leading.indicators.incurred.lr[,"Variable"]]
write.csv(lead.ind.inc.lr.matrix, "lead.ind.inc.lr.matrix.csv")
#write.csv(x.scaled, "x.csv")

#identify the longest lag time to adjust dependent (y) variable
longest.lag <- max(abs(as.numeric(leading.indicators.incurred.lr[,"Lag"])))
shortest.lag <- min(abs(as.numeric(leading.indicators.incurred.lr[,"Lag"])))

#save variable names for the significant leading indicators
leading.ind.headings <- colnames(lead.ind.inc.lr.matrix)

##adjust the matrix for each variable's corresponding lag time
#lead.ind.inc.lr.matrix[,1]
lead.ind.inc.lr.matrix <- apply(lead.ind.inc.lr.matrix, 2, rev) #reverse order of variables
#lead.ind.inc.lr.matrix[,1] #check to be sure reverse ording worked
temp.lag.adj <- NULL
lead.ind.lag.adjusted <- NULL
#i=1
for(i in 1:length(leading.ind.headings)){
#lead.ind.inc.lr.matrix[,i]
  temp.lag.adj <- lead(lead.ind.inc.lr.matrix[,i], abs(as.numeric(leading.indicators.incurred.lr[i,"Lag"])))
  lead.ind.lag.adjusted <- cbind(lead.ind.lag.adjusted, temp.lag.adj)
}

#lead.ind.lag.adjusted[,1]

#reverse back the order of the dependent variable training set, after it was adjusted for lags
lead.ind.lag.adjusted <- apply(lead.ind.lag.adjusted, 2, rev) #reverse order of variables

#lead.ind.lag.adjusted[,1]
#length(lead.ind.lag.adjusted[1,])

#add back the column headers
tidy.colnames <- make.names(leading.ind.headings, unique=TRUE)
colnames(lead.ind.lag.adjusted) <- tidy.colnames
write.csv(lead.ind.lag.adjusted, "lead.ind.lag.adjusted.including.NAs.csv")

y.leading.nonwhitened.unscaled <- drop_na(data.frame(lead(y.unscaled, longest.lag)))
rownames(y.leading.nonwhitened.unscaled) <- NULL
y.leading.nonwhitened.unscaled <- as.matrix(y.leading.nonwhitened.unscaled)
colnames(y.leading.nonwhitened.unscaled) <- make.names(colnames(y.leading.nonwhitened.unscaled))


#remove the NA's, which also has the effect of adjusting the number of rows to be less than the greatest lag
lead.ind.lag.adjusted <- drop_na(data.frame(lead.ind.lag.adjusted))

x.leading.whitened.scaled <- data.frame(lead.ind.lag.adjusted)
rownames(x.leading.whitened.scaled) <- NULL #remove rownames, which get messed up when using rev
x.leading.whitened.scaled <- as.matrix(x.leading.whitened.scaled)
write.csv(x.leading.whitened.scaled, "x.leading.whitened.scaled.csv")

#append dependent variable: only used for modeling algorithms that need y ~ x format, instead of x, y
lead.ind.lag.adjusted <- data.frame(x.leading.whitened.scaled, y.leading.nonwhitened.unscaled)
write.csv(lead.ind.lag.adjusted, "lead.ind.lag.adjusted.csv")

### generate model ###

#model.lm <- lm(y ~ ., data=lead.ind.lag.adjusted)
#model.glm <- glm(y ~ ., data=lead.ind.lag.adjusted)
#cvfitlm <- cv.glmnet(x = lead.ind.lag.adjusted[,-length(lead.ind.lag.adjusted)],  y = as.numeric(lead.ind.lag.adjusted[,"y"]), family = "gaussian", alpha = 0.5)

##fit an ordinal lasso model
#?glmnetcr
#cor(x.leading.whitened.scaled, y.leading.nonwhitened.unscaled)
##ordinal glmnet
#glmnet.fit <- glmnetcr(x.leading.whitened.scaled, y.leading.nonwhitened.unscaled, nlambda=10000000, alpha=0)
#BIC.model <- select.glmnetcr(glmnet.fit)
#estimates <- coef(glmnet.fit, s = BIC.model)
#fitted(glmnet.fit, s = select.glmnetcr(cvfitlm))
#plot(glmnet.fit ,xvar="step", type="bic")

##basic glmnet
#glmnet.fit <- glmnet(x.leading.whitened.scaled, y.leading.nonwhitened.unscaled, nlambda=10000000, alpha=1, lambda.min.ratio=0.5)
#str(glmnet.fit)
#plot(glmnet.fit)

##cross validated glmnet
cv.glmnet.fit <- cv.glmnet(x.leading.whitened.scaled, y.leading.nonwhitened.unscaled, family = "gaussian", alpha = 1)
str(cv.glmnet.fit)

pdf("GLMNET_MSE.pdf")
plot(cv.glmnet.fit)
dev.off()

@

\section{GLM Net Model}

Figure~\ref{fig:glmnetmse} shows...

\begin{figure}[h!]
\begin{center}
\includegraphics[width=\linewidth]{GLMNET_MSE.pdf}
\caption{Mean Squared Error for different values of lambda}
\label{fig:glmnetmse}
\end{center}
\end{figure}

<<echo=FALSE, results=hide>>=

#get coefficients from GLMNET
cvfitlm.coef.lambda.1se.prewhitened <- coef(cv.glmnet.fit, s = "lambda.1se")

#send coefficient list to a text file
capture.output(print(cvfitlm.coef.lambda.1se.prewhitened), file="GLMNET_Coefficients.txt")

@

\lstinputlisting{GLMNET_Coefficients.txt}

<<echo=FALSE, results=hide>>=

#str(cvfitlm.coef.lambda.1se.prewhitened)
#model.names <- as.data.frame(cvfitlm.coef.lambda.1se.prewhitened@Dimnames[1])
#variable.indices <- as.data.frame(cvfitlm.coef.lambda.1se.prewhitened@i)

#str(model.names)
#str(variable.indices)
#model.names[1]
#model.names[variable.indices]

#model.coefs <- cvfitlm.coef.lambda.1se.prewhitened@x
#str(model.names)
#length(model.names[1])
#length(model.coefs)

##use the model results from the pre-whitened data to select the variables to use in the forecast
#coef.indices.lambda.1se.prewhitened <- which(cvfitlm.coef.lambda.1se.prewhitened@x != 0)
#selected.indices <- coef.indices.lambda.1se.prewhitened[-1] #-1 is to omit the y-intercept variable
#x.selected <- x_new[,selected.indices]
##x_new_colnames[selected.indices] #as a test

##create forecast set matrix based on each variable's corresponding lag time
#lead.ind.inc.lr.matrix[,2]
forecast.raw.matrix <- apply(lead.ind.inc.lr.matrix, 2, rev) #reverse order of variables
#forecast.raw.matrix[,1] #check to be sure reverse ording worked
#lead.ind.inc.lr.matrix[,1] #compare
temp.lag.adj <- NULL
forecast.set <- NULL
#i=1
for(i in 1:length(leading.ind.headings)){
#lead.ind.inc.lr.matrix[,i]
  temp.lag.adj <- lead(lead.ind.inc.lr.matrix[,i], abs(as.numeric(leading.indicators.incurred.lr[i,"Lag"])+shortest.lag))
  forecast.set <- cbind(forecast.set, temp.lag.adj)
}

#forecast.set[,1]

#reverse back the order of the dependent variable training set, after it was adjusted for lags
forecast.set <- apply(forecast.set, 2, rev) #reverse order of variables
#forecast.set[,1]

#add back the column headers
tidy.colnames <- make.names(leading.ind.headings, unique=TRUE)
colnames(forecast.set) <- tidy.colnames
write.csv(forecast.set, "forecast_set_NAs.csv")

#remove the NA's, which also has the effect of adjusting the number of rows to be less than the greatest lag
forecast.set <- drop_na(data.frame(forecast.set))
forecast.matrix <- as.matrix(forecast.set)
#forecast.set[,1]
write.csv(forecast.set, "forecast_set.csv")

#x.leading.whitened.scaled <- data.frame(lead.ind.lag.adjusted)
#rownames(x.leading.whitened.scaled) <- NULL #remove rownames, which get messed up when using rev
#x.leading.whitened.scaled <- as.matrix(x.leading.whitened.scaled)
#write.csv(x.leading.whitened.scaled, "x.leading.whitened.scaled.csv")



##make predictions
predictions <- predict(cv.glmnet.fit, newx=forecast.matrix, type="response", s = "lambda.1se")
?predict.cv.glmnet

##create plot of the glmnet forecast
y.actual <- y.leading.nonwhitened.unscaled
forecast.ahead <- as.matrix(forecast.ahead)
forecast.ahead <- drop_na(data.frame(predictions[length(y.actual)+1:length(predictions)]))

pdf("glmnet_forecast.pdf")
plot(c(y.actual, forecast.ahead), type='p', main = "Forecast", ylab="Values", xlab="Time")
lines(y.actual, col = "black", lty = 1, lwd = 2)
lines(predictions, col = "green", lty = 1, lwd = 2)
legend("topleft", c("Forecast", "Actual"), col = c("green", "black"), text.col = "black", lty = c(1, 1, 1, 1, 2), lwd = c(2, 2, 2, 1, 2), merge = TRUE, bg = 'gray90', cex = .75)
dev.off()

@

\section{GLM Net Model}

Figure~\ref{fig:glmnetforecast} shows...

\begin{figure}[h!]
\begin{center}
\includegraphics[width=\linewidth]{glmnet_forecast.pdf}
\caption{Forecast versus actual, based on cross validated GLM Net}
\label{fig:glmnetforecast}
\end{center}
\end{figure}


<<echo=FALSE, results=hide>>=

colnames(forecast.ahead) <- "Forecast (ascending order if more than one value)"
kable(forecast.ahead)
write.csv(file="y-prediction.csv", x=predict.ahead)


###############################################################
##Forcast with BVAR method

Y <- lead.ind.lag.adjusted
#colnames(Y) <- col.names
nms <- dependent.variable
Y <- as.matrix(Y)

# Fit a Basic VAR-L(3,4) on simulated data
T1=floor(nrow(Y)/3)
T2=floor(2*nrow(Y)/3)
#?constructModel
#m1=constructModel(Y,p=4,struct="Basic",gran=c(20,10),verbose=FALSE,IC=FALSE,T1=T1,T2=T2,ONESE=TRUE)
#m1=constructModel(Y,p=4,struct="Tapered",gran=c(50,10),verbose=FALSE,T1=T1,T2=T2,IC=FALSE)
#plot(m1)
#results=cv.BigVAR(m1)
#plot(results)
#predict(results,n.ahead=1)

#SparsityPlot.BigVAR.results(results)

#str(results)
#results@preds
#results@alpha
#results@Granularity
#results@Structure
#results@lagmax
#results@Data
#plot(results@Data)

#install.packages("devtools")
#library(devtools)
#install_github("gabrielrvsc/HDeconometrics")

###################################
#The above, BigVAR package will not handle data sets this wide.  Trying the
#Bayesian Vector Auto Regression (BVAR) algorithm

###################################
##Perform analysis on pre-whitened data

# = load package and data = #
#install.packages("HDeconometrics")
#library(HDeconometrics)
#data("voldata")

# = Break data into in and out of sample to test model accuracy= #
#Yin=voldata[1:5499,]
#Yout=voldata[-c(1:5499),]
Yin = Y[1:T2,]
Yout = Y[(T2+1):(T1+T2),]

# = Run models = #
# = OLS = #
#modelols=HDvar(Yin,p=2) # takes a while to run
#predols=predict(modelols,h=2)

# = BVAR = #
#?lbvar
#?predict
## == This example uses the Brazilian inflation data from
#Garcia, Medeiros and Vasconcelos (2017) == ##
#data("BRinf")
#Y=BRinf[,1:59]# remove expectation variables
#modelB=lbvar(Y,p=4)

# take a look at the coefficients
#eq=coef(modelB,type="equation")
#block=coef(modelB,type="block")
#block$Lag1

modelbvar=lbvar(Yin, p = 2)

# take a look at the coefficients
eq=coef(modelbvar,type="equation")
block=coef(modelbvar,type="block")
block$Lag1

predbvar=predict(modelbvar,h=2)

# = 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])),predols[,k]))
lines(c(rep(NA,length(Y[,k])),predbvar[,k]))
abline(v=length(Y[,k]),lty=2,col=4)
#legend("topleft",legend="BVAR",col=2,lty=1,lwd=1,seg.len=1,cex=1,bty="n")
dev.off()

# = Overall percentual error = #
#MAPEols=abs((Yout-predols)/Yout)*100
#MAPEbvar=abs((Yout-predbvar)/Yout)*100
#matplot(MAPEols,type="l",ylim=c(0,80),main="Overall % error",col="lightsalmon",ylab="Error %")
#aux=apply(MAPEbvar,2,lines,col="lightskyblue1")
#lines(rowMeans(MAPEols),lwd=3,col=2,type="b")
#lines(rowMeans(MAPEbvar),lwd=3,col=4,type="b")
#legend("topleft",legend=c("OLS","BVAR"),col=c(2,4),lty=1,lwd=1,seg.len=1,cex=1,bty="n")

# = Influences = #
#aux=modelbvar$coef.by.block[2:23]
#impacts=abs(Reduce("+", aux ))
#diag(impacts)=0
#I=colSums(impacts)
#R=rowSums(impacts)
#par(mfrow=c(2,1))
#barplot(I,col=rainbow(30),cex.names = 0.3, main = "Most Influent")
#barplot(R,col=rainbow(30),cex.names = 0.3, main = "Most Influenced")

#pdf(file=paste("plots/", dependent.variable, "_barchart.pdf",sep=""))
#aux=modelbvar$coef.by.block
#impacts=abs(Reduce("+", aux ))
#diag(impacts)=0
#I=colSums(impacts)
#R=rowSums(impacts)
#par(mfrow=c(2,1))
#barplot(I,col=rainbow(30),cex.names = 0.3, main = "Most Influent")
#barplot(R,col=rainbow(30),cex.names = 0.3, main = "Most Influenced")
#dev.off()

###################################
##Perform analysis on NON pre-whitened data

Y <- cbind.data.frame(x.scaled, y.unscaled)
#colnames(Y) <- col.names
head(Y)
nms <- dependent.variable
Y <- as.matrix(Y)

# Fit a Basic VAR-L(3,4) on simulated data
T1=floor(nrow(Y)/3)
T2=floor(2*nrow(Y)/3)



# = Break data into in and out of sample to test model accuracy= #
Yin = Y[1:T2,]
Yout = Y[(T2+1):(T1+T2),]

# = Run models = #
# = OLS = #
#modelols=HDvar(Yin,p=2) # takes a while to run
#predols=predict(modelols,h=2)

# = BVAR = #
#?lbvar
#?predict
modelbvar=lbvar(Yin, p = 2, delta = 0)
predbvar=predict(modelbvar,h=2)

# = Forecasts of the volatility = #
k=paste(dependent.variable)
pdf(file=paste("plots/", dependent.variable, "_forecast_not_whitened.pdf",sep=""))
plot(c(Y[,k],predbvar[,k]),type="l", main=paste(dependent.variable, "Not Whitened"))
#lines(c(rep(NA,length(Y[,k])),predols[,k]))
lines(c(rep(NA,length(Y[,k])),predbvar[,k]))
abline(v=length(Y[,k]),lty=2,col=4)
#legend("topleft",legend="BVAR",col=2,lty=1,lwd=1,seg.len=1,cex=1,bty="n")
dev.off()

# = Overall percentual error = #
#MAPEols=abs((Yout-predols)/Yout)*100
#MAPEbvar=abs((Yout-predbvar)/Yout)*100
#matplot(MAPEols,type="l",ylim=c(0,80),main="Overall % error",col="lightsalmon",ylab="Error %")
#aux=apply(MAPEbvar,2,lines,col="lightskyblue1")
#lines(rowMeans(MAPEols),lwd=3,col=2,type="b")
#lines(rowMeans(MAPEbvar),lwd=3,col=4,type="b")
#legend("topleft",legend=c("OLS","BVAR"),col=c(2,4),lty=1,lwd=1,seg.len=1,cex=1,bty="n")

#bar chart
#pdf(file=paste("plots/", dependent.variable, "_barchart_not_whitened.pdf",sep=""))
#aux=modelbvar$coef.by.block
#impacts=abs(Reduce("+", aux ))
#diag(impacts)=0
#I=colSums(impacts)
#R=rowSums(impacts)
#par(mfrow=c(2,1))
#barplot(I,col="blue",cex.names = 0.3, main = "Most Influential")
#barplot(R,col="blue",cex.names = 0.3, main = "Most Influenced")
#dev.off()

###############################################################