CPI_Sweave_Template.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{Consumer Price Index 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 portend a possible increase in inflation.  For example, the U.S. national debt to GDP has rarely been higher, and has caused great concern for many.  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 inflation, as well.  This analysis examines these variables, and their combined relationship to inflation, as measured by CPI.  The data in this study goes back to 1961, and extends through 2019.

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.  As older, retired individuals tend to purchase fewer consumer goods, this has helped inflation to remain very low.  

The difference between the number of people entering age 65, the standard retirement age, compared to those entering age 45, the age many reach their peak earnings growth potential, can provide a useful indicator as to the net impact of age demographics on inflation.  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())
#Sys.which("pdflatex")
#install.packages("orderedLasso")

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("C:\\Users\\P2889196\\Documents\\Code\\mygithub\\BayesianVectorAutoRegression")
##########################

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

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

#colnames(data) <- c("CPI","shortTermDebtToLongTerm","debtToGDP","GDP")

##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:59,3:5]
births.2009 <- apply(births.trunc, 2, rev)
head(births.2009)

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

col.names <- colnames(SeriesData)

x <- SeriesData[, !(names(SeriesData) %in% dependent.variable)]
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 - not used
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
#i=1
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.unscaled)
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 the Consumer Price Index}

Before examining the combined, multivariate impact of all the variables on changes in the Consumer Price Index, this section will perform a univariate analysis to create a forecast on nothing more than the CPI 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:ArimaCPI} shows that this method estimates that the year to year trend in the CPI will be declining for the next five years, beginning in 2018.  

As an ARIMA forecast is univariate, this method does not take into account any of the previously mentioned economic and demographic variables.  It only assumes that past patterns in the trend of the CPI will continue.  In the next section, we will evaluate the potential impact of the previously mentioned economic and demographic variables.

\begin{figure}[!ht]
\begin{center}
\includegraphics[width=\linewidth]{plots/CPI_arima.pdf}
\caption{Univariate ARIMA forecast of the Consumer Price Index}
\label{fig:ArimaCPI}
\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.unscaled)
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+1),]

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

# = BVAR Forecasts = #
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 BVAR}

In this section, we produce a multivariate forecast of the CPI, 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:BVARForecastCPI} shows a forecast of the CPI using the BVAR method.  The black trend line to the right of the blue, dashed, vertical line, indicates the forecasted trend in the CPI for five years, after 2017.  The forecast suggests that annual CPI will largely increase on the whole.

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

\section{Influencers}

Figure~\ref{fig:BVARInfluencersCPI} 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 CPI 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:BVARInfluencersCPI} that all the predictor variables tend to be relatively strong influencers of the velocity of M2 money supply, as shown in the variable, m2\_velocity.  (M2 money supply includes the cash and deposits represented by M1 money supply, plus what is known as near money, which includes such things as savings and money market accounts, and other cash deposits.)  However, the predictor variables are not as strong as influencers of the Consumer Price Index.  There may be other variables that are more predictive, and could increase forecast confidence.

It might be surprising to some, that the debt to GDP ratio does not appear as a more prominent influencer, given all the media attention the national debt has received.  However, there have been other years in history that saw worse debt to GDP ratios, and it did not necessarily drive high inflation.  Also, while some have suggested that the millenial generation may influence inflation upwards as they begin to set up households, and increase demand for all the goods and servies that go with that endeavor, it does not appear that this is a primary influencer of inflation this time around.  It may be that the Baby Boomer generation is still a primary counter-acting factor.

Since the velocity of M2 money supply is so affected by the other predictor variables, it might be interesting to see how a forecast of M2 might be affected.  It is easy to imagine how inflation, demographics, and debt factors would affect the velocity of money moving in and out of the highly liquid accounts represented by M2 money supply.  It is not certain whether the direction of velocity would move upward or downward in such a forecast, however. 

<<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}[!ht]
\begin{center}
\includegraphics[width=\linewidth]{plots/CPI_barchart_not_whitened.pdf}
\caption{Main inflencers of the Consumer Price Index - Bayesian Vector Auto Regression (BVAR)}
\label{fig:BVARInfluencersCPI}
\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 inflation.  Many have predicted an increase in inflation, 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, do indeed point to the possibility of increasing inflation.

\end{document}