DJI_Sweave_Template copy.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{Dow Jones 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{Dow Jones Forecast Variables and Methods}

Following is a brief analysis, to forecast the Consumer Price Index, as a measure of inflation.  It is not meant to be a comprehensive, broad treatment of the many variables that might affect inflation.  Rather, the analysis focuses on some economic trends that have been in the news lately, and that are sometimes said to portent a possible crash in the markets.  For example, the U.S. national debt to GDP has rarely been higher, and has caused great concern for many investors.  Also, the rate of inflation has just begun to increase, after years of being at record lows.  Lastly among economic variables, the velocity of money is sometimes said to be related to market performance.  This analysis examines these variables, and their combined relationship to the year to year changes in the Dow Jones Industrial Index.  The data in this study goes back to 1961, and extends through 2017.

In addition to the economic variables, demographic variables were also considered.  As a proxy for age distributions, the number of births by year, and also the birth rate by year, were included along with the previously mentioned economic variables.  It is well known that the population of the U.S. has been aging.  This trend has contributed to rising demand for fixed income investments, which has driven interest rates to very low levels.  Also, as older, retired individuals tend to purchase fewer consumer goods, inflation for these items as also been very low.  This has an impact toward lower velocity of money as well.  As retirees have been known to rebalance investment portfolios away from equities, and toward fixed income investments, large numbers of retirees have been seen to have a negative impact on the markets.  

Conversely, when there are large numbers of people in their peak earnings years, which is approximately age 45 on average, it tends to have a positive impact on the markets.  We saw such a period in the 1990s, when a large proportion of Baby Boomers were in their mid forties.  

The difference between the number of people in these two groups, those entering age 65 compared to those entering age 45, can provide a useful indicator as to the net impact of age demographics on the markets.  This variable was also included in this analysis, to determine how demographic forces interact with the other key economic variables in this forecast.


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

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

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

##########################
#Type 'TRUE' if you want to include an offset in the GLM calculation, FALSE otherwise.
include.offset <- FALSE
##########################

##########################
#Type the column name of the variable you would like to use as an offset, if any.
#offset.variable <- "UnitCount"
##########################

##########################
#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="annual", start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)

#Nonfinancial corporate business; short-term debt as a percentage of total debt, Annual
shortTermDebtToLongTerm <- Quandl("FED/FL104140006_A", api_key="DJGcfzQc5RYP1JSycMBv", collapse="annual", 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="annual", start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)


GDP <- Quandl("FED/FU086902001_A", api_key="DJGcfzQc5RYP1JSycMBv", collapse="annual", 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="annual", start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)


m2Velocity <- Quandl("FRED/M2V", api_key="DJGcfzQc5RYP1JSycMBv", collapse="annual", start_date="1960-12-31", end_date="2017-12-31", type="raw", order="asc", force_irregular=TRUE)

raw.ts <- cbind(m1Velocity, m2Velocity, GDP[-1,], debtToGDP[-1,], shortTermDebtToLongTerm[-1,], CPI[-1,])

raw.ts <- cbind(m1Velocity, m2Velocity, GDP[-1,], debtToGDP[-1,], shortTermDebtToLongTerm[-1,], CPI[-1,])

#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","shortTermDebtToLongTerm","debtToGDP","GDP","m1Velocity","m2Velocity")

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

SeriesData <- cbind(births.2009, data, DJIChg)
SeriesData <- SeriesData[,-1] #remove the "Year" column
head(SeriesData)

col.names <- colnames(SeriesData)

x <- SeriesData[, !(names(SeriesData) %in% dependent.variable)]
head(x)

#scale the dependent variable
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 <- colnames(x)


## ----ARIMA, results='asis'-----------------------------------------------
#i=20
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])
  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)

@

\section{Arima Univariate Forecast of Dow Jones Industrial Average}

Before examining the combined, multivariate impact of all the variables on changes in the Dow Jones Industrial Average, this section will perform a univariate analysis to create a forecast on nothing more than the DJIA data alone.  The method utilized is an auto-ARIMA, or automatic Auto Regressive Integrated Moving Average.  This brief document will not go into the details as to how the ARIMA method works.  However, perhaps it will suffice to say that it examines the following.

\begin{itemize}
\item Auto Regressive trends - the AR in ARIMA: measures a trends correlation to itself, such as with seasonality, or in the case of annual data, possible cyclicality
\item Integration - the I in ARIMA: measures the degree to which subracting each preceding time period's value removes the overall trend, also known as detrending.
\item Moving Average - the MA in ARIMA: measures the extent to which a trend can be explained by a moving average pattern
\end{itemize}

The auto.arima function in R estimates each of these ARIMA parameters, and uses them to calculate a forecast.  The plot in Figure~\ref{fig:ArimaDJI} shows that this method estimates that the year to year change in the DJIA will be flat for the next five years, beginning in 2018.  Also, at the top of the chart, we see that the ARIMA parameters chosen by the algorithm are "(0,0,0)".  This indicates no discernable AR, I, or MA effects. 

As an ARIMA forecast is univariate, this method does not take into account any of the previously mentioned economic and demographic variables.  The ARIMA analysis was performed here to see whether there were any affects within the DJIA trend itself that should be taken into account, which might significantly affect the forecast.  There do not appear to be any strong ARIMA affects that would affect our forecast either up, or down.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=\linewidth]{plots/DJIChg_arima.pdf}
\caption{Univariate ARIMA forecast of the Dow Jones Industrial Index}
\label{fig:ArimaDJI}
\end{center}
\end{figure}

<<echo=False, results=hide>>=
## ------------------------------------------------------------------------
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()
}


## ------------------------------------------------------------------------
#i=2
##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
#i=1
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()
  }
}


## ----analytical_dataset, results='hide'----------------------------------

#x
#y
#x.arima.residuals <- apply(x.arima.residuals, 2, rev)
#y.arima.residuals <- rev(y.arima.residuals)
#time.increments <- rev(time.increments)

#x2 <- subset(x, select= -c(X.Dividends, X.From.the.rest.of.the.world.6.))


Y <- cbind.data.frame(x, y)
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
modelbvar=lbvar(Yin, p = 5, delta = 0.5)
predbvar=predict(modelbvar,h=5)

# = 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 = 5, delta = 0.5)
predbvar=predict(modelbvar,h=5)

# = 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")

@

\section{Multivariate Forecast with Bayesian Vector Auto Regression - BVAR}

In this section, we produce a multivariate forecast of the DJIA, utilizing the previously mentioned economic and demographic variables.  Older methods of performing multivariate time series analysis involved searching for top correlations between the predictor variables, and the target variable we are trying to predict.  It largely ignored the impact of predictor variables on each other.  A newer method for accounting for these cross variable impacts, is called Bayesian Vector Auto Regression, or BVAR.

Figure~\ref{fig:BVARForecast} shows a forecast of the DJIA using the BVAR method.  The black trend line to the right of the blue, dashed, vertical line, indicates the forecasted trend in DJIA annual percentage growth for five years, after 2017.  The forecast suggests that annual growth will moderate, diminishing year over year.  The growth rate, however, is not predicted to be negative on a year over year basis.  Crashes in the short run as always possible.  However, these economic indicators, and demographic age factors taken together, do not indicate strong pressures for a market crash over the next five years, despite much media hyperbole to the contrary.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=\linewidth]{plots/DJIChg_forecast.pdf}
\caption{Forecast of the Dow Jones Industrial Index, utilizing predictor variables that are scaled, and pre-whitened - Bayesian Vector Auto Regression (BVAR)}
\label{fig:BVARForecast}
\end{center}
\end{figure}

\section{Influencers}

Figure~\ref{fig:BVARInfluencers} shows the top variables in the BVAR forecast model, as to their influence on each other.  Recall that the BVAR method takes into account not only each variable's impact on the target variable, the DJIA in this case, but also on each other.  The BVAR also examines the relationships between the variables, adjusting for many possible lags of the data over time.  This enables finding the variables that tend to be leading indicators.  This is a form of causality analysis, which investigates whether the strongest relationship, akin to the concept of correlation, occurs when the movement of one variable, precedes a corresponding movement on the part of another.

We can see in Figure~\ref{fig:BVARInfluencers} that all the predictor variables not only tend to be relatively strong influencers of each other, but especially of the DJIA.  The DJIA, on the other hand, is not a strong influencer of the other variales.  

This analysis also shows that the trend in the change in the DJIA is not a strong influencer of itself.  In other words, the past trend of the year over year trend in the DJIA is not a good predictor of the future.  We also saw this in the univariate ARIMA analysis as well.

<<name="influencers", echo=False, fig=False>>=
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()

@

\begin{figure}[h!]
\begin{center}
\includegraphics[width=\linewidth]{plots/DJIChg_barchart_not_whitened.pdf}
\caption{Main inflencers of the Dow Jones Industrial Index - Bayesian Vector Auto Regression (BVAR)}
\label{fig:BVARInfluencers}
\end{center}
\end{figure}

\section{Summary}

This analysis examined the impact of some economic and demographic trends that have gotten attention in the media, as to their influence on the stock market.  Many have predicted a collapse in the stock market, as a result of some of these trends.  While no analysis can guarantee against such a collapse, and there are many other variables involved with market performance than these few, this analysis suggests that recent trends in these variables, when taken together, have not reached levels that would make a severe decline in the markets inevitable.  However, the resulting forecast from this mulivariate analysis does suggest that the recent double digit growth rates of the market are not likely to be sustainable.  While annual growth rates are likely to remain positive, they will likely diminish year over year for at least the next five years.

\end{document}