Outline

Probability Models for Cycles

Autoregression Model

Properties of Autoregression Models

Lag Polynomials

Stationarity Conditions

Integrated Processes

Testing Approach

The Stakes: The “Unit Root Wars”

Application: Macroeconomic Forecasts

Stationarity and AR order choice results

#Libraries
library(fredr) # Data from FRED API
library(fpp2) #Forecasting and Plotting tools
library(vars) #Vector Autoregressions
library(knitr) #Use knitr to make tables
library(kableExtra) #Extra options for tables
library(dplyr) #Data Manipulation
library(tseries) #Time series functions including stationarity tests
library(gridExtra)  #Graph Display

# Package "BMR" for BVAR estimation is not on CRAN, but is instead maintained by an individual 
# It must be installed directly from the Github repo: uncomment the following code to do so

# library(devtools) #Library to allow downloading packages from Github
# install_github("kthohr/BMR")
# library(BMR) #Bayesian Macroeconometrics in R

# Note that if running this code on Kaggle, internet access must be enabled to download and install the package
# If installed locally, there may be difficulties due to differences in your local environment (in particular, versions of C++)
# For this reason, relying local installation is not recommended unless you have a spare afternoon to dig through help files

#An alternative, similar library, is BVAR: it is on CRAN
library(BVAR) #Bayesian Vector Autoregressions

##Obtain and transform NIPA Data (cf Lecture 08)

fredr_set_key("8782f247febb41f291821950cf9118b6") #Key I obtained for this class

## Load Series: Series choices and names as in Litterman (1986)

RGNP<-fredr(series_id = "GNPC96",
           observation_start = as.Date("1971-04-01"),
           observation_end = as.Date("2020-01-01"),
           vintage_dates = as.Date("2021-03-29"),
           units="cch") #Real Gross National Product, log change

INFLA<-fredr(series_id = "GNPDEF",
           observation_start = as.Date("1971-04-01"),
           observation_end = as.Date("2020-01-01"),
           vintage_dates = as.Date("2021-03-29"),
           units="cch") #GNP Deflator, log change

UNEMP<-fredr(series_id = "UNRATE",
           observation_start = as.Date("1971-04-01"),
           observation_end = as.Date("2020-01-01"),
           vintage_dates = as.Date("2021-03-29"),
           frequency="q") #Unemployment Rate, quarterly

M1<-fredr(series_id = "M1SL",
           observation_start = as.Date("1971-04-01"),
           observation_end = as.Date("2020-01-01"),
           vintage_dates = as.Date("2021-03-29"),
           frequency="q",
           units="log") #Log M1 Money Stock, quarterly

INVEST<-fredr(series_id = "GPDI",
           observation_start = as.Date("1971-04-01"),
           observation_end = as.Date("2020-01-01"),
           vintage_dates = as.Date("2021-03-29"),
           units="log") #Log Gross Domestic Private Investment

# The 4-6 month commercial paper rate series used in Litterman (1986) has been discontinued: 
# For sample continuity, we merge the series for 3 month commercial paper rates from 1971-1997 with the 3 month non-financial commercial paper rate series
# This series also has last start date, so it dictates start date for series

CPRATE1<-fredr(series_id = "WCP3M",
           observation_start = as.Date("1971-04-01"),
           observation_end = as.Date("1996-10-01"),
           vintage_dates = as.Date("2021-03-29"),
           frequency="q") #3 Month commercial paper rate, quarterly, 1971-1997

CPRATE2<-fredr(series_id = "CPN3M",
           observation_start = as.Date("1997-01-01"),
           observation_end = as.Date("2020-01-01"),
           vintage_dates = as.Date("2021-03-29"),
           frequency="q") #3 Month AA nonfinancial commercial paper rate, quarterly, 1997-2018

CPRATE<-full_join(CPRATE1,CPRATE2) #Merge 2 series to create continuous 3 month commercial paper rate series from 1971-2018

CBI<-fredr(series_id = "CBI",
           observation_start = as.Date("1971-04-01"),
           observation_end = as.Date("2020-01-01"),
           vintage_dates = as.Date("2021-03-29")) #Change in Private Inventories

#Format the series as quarterly time series objects, starting at the first date
rgnp<-ts(RGNP$value,frequency = 4,start=c(1971,2),names="Real Gross National Product") 
infla<-ts(INFLA$value,frequency = 4,start=c(1971,2),names="Inflation")
unemp<-ts(UNEMP$value,frequency = 4,start=c(1971,2),names="Unemployment")
m1<-ts(M1$value,frequency = 4,start=c(1971,2),names="Money Stock")
invest<-ts(INVEST$value,frequency = 4,start=c(1971,2),names="Private Investment")
cprate<-ts(CPRATE$value,frequency = 4,start=c(1971,2),names="Commercial Paper Rate")
cbi<-ts(CBI$value,frequency = 4,start=c(1971,2),names="Change in Inventories")


#Express as a data frame
macrodata<-data.frame(rgnp,infla,unemp,m1,invest,cprate,cbi)

nlags<-6 # Number of lags to use
nseries<-length(macrodata[1,]) #Number of series used
#Evaluate Stationarity of series by Phillips-Perron Tests, as well as KPSS tests
# accounting for constants and possible linear trend
#Use auto.arima to choose AR order after KPSS test without trend
stationaritytests<-list()
PPtests<-list()
KPSStests<-list()
ARIstatmodels<-list()
Integrationorder<-list()
ARorder<-list()
for (i in 1:nseries){
  PPtests[i]<-PP.test(macrodata[,i])$p.value
  KPSStests[i]<-kpss.test(macrodata[,i],null="Trend")$p.value
  ARIstatmodels[[i]]<-auto.arima(macrodata[,i],max.q=0,seasonal=FALSE) #Apply auto.arima set to (nonseasonal) ARI only
  Integrationorder[i]<-ARIstatmodels[[i]]$arma[6] #Integration order chosen (uses KPSS Test)
  ARorder[i]<-ARIstatmodels[[i]]$arma[1] #AR order chosen (uses AIC)
}
Series<-c("Real GNP Growth","Inflation","Unemployment","Money Stock","Private Investment","Commercial Paper Rate","Change in Inventories")
stationaritytests<-data.frame(Series,as.numeric(PPtests),as.numeric(KPSStests),as.numeric(Integrationorder),as.numeric(ARorder))

rownames(stationaritytests)<-Series

colnames(stationaritytests)<-c("Series","PP p-value","KPSS p-value","Integration Order","AR Order")

stationaritytests %>%
  mutate(
    Series = row.names(.),
    `PP p-value` = cell_spec(`PP p-value`,color = ifelse(`PP p-value`<0.05, "black", "red")),
    `KPSS p-value` = cell_spec(`KPSS p-value`,color = ifelse(`KPSS p-value`<0.05, "red", "black")),
    `Integration Order`= cell_spec(`Integration Order`,color = ifelse(`Integration Order`==1, "red", "black")),
    `AR Order`=`AR Order`
  ) %>%
kable(escape=F,
  caption="Chosen Models and P-Values for Stationarity Tests: Red indicates Unit Root") %>%
  kable_styling(bootstrap_options = "striped")
Chosen Models and P-Values for Stationarity Tests: Red indicates Unit Root
Series PP p-value KPSS p-value Integration Order AR Order
Real GNP Growth Real GNP Growth 0.01 0.1 0 2
Inflation Inflation 0.01 0.01 1 4
Unemployment Unemployment 0.395363350741117 0.0512422827489796 1 1
Money Stock Money Stock 0.918305692612384 0.01 1 2
Private Investment Private Investment 0.470464420263501 0.01 1 1
Commercial Paper Rate Commercial Paper Rate 0.0771149514500132 0.0180787556682605 1 3
Change in Inventories Change in Inventories 0.01 0.1 1 0

Real GNP Growth: Series, ACF, PACF, Estimated AR(2) Roots

rgnpplots<-list()

rgnpplots[[1]]<-autoplot(rgnp)+labs(x="Date",y="Percent Growth",title="Real GNP Growth")
rgnpplots[[2]]<-ggAcf(rgnp)+labs(title="Autocorrelation Function")
rgnpplots[[3]]<-ggPacf(rgnp)+labs(title="Partial Autocorrelation Function")
rgnpplots[[4]]<-autoplot(ARIstatmodels[[1]])
grid.arrange(grobs=rgnpplots,nrow=2,ncol=2)

Priors for Autoregression Models

Variance Priors

Minnesota Priors

The Multivariate Case: Vector Autoregression Model

Multivariate VAR Priors

Application: Litterman (1986) Macroeconomic Forecasts

Forecasts from ARI (Blue) and BVAR (Red)

#Label observations with dates
rownames(macrodata)<-M1$date

#Set priors: use hierarchical approach, with default scales, to account for uncertainty
#Note that defaults are for series in levels, while some here are growth rates
mn <- bv_minnesota(lambda = bv_lambda(mode = 0.2, sd = 0.4, min = 0.0001, max = 5),
                   alpha = bv_alpha(mode = 2), 
                   var = 1e07)

#Ghost samples to set priors on initial conditions: 
# SOC shifts towards independent unit roots
# SUR shifts toward shared unit root
# Hierarchical approach allows strength to be determined by data
#soc <- bv_soc(mode = 1, sd = 1, min = 1e-04, max = 50)
#sur <- bv_sur(mode = 1, sd = 1, min = 1e-04, max = 50)

priors <- bv_priors(hyper = c("lambda"), #Choose lambda by hierarchical method, leave others fixed at preset
                    mn = mn) #Set Minnesota prior as defined above

#Add options soc = soc, sur = sur to set initial observations priors
# For this data, hierarchical estimation sets fairly smal values for these, and resulting forecasts are about the same 
# except with slightly narrower intervals for far future forecasts

#Set up MCMC parameters: use default initialization except that scaling is adjusted automatically in burnin
#If this is not enough, may need to adjust manually: check diagnostics
mh <- bv_metropolis(adjust_acc = TRUE, acc_lower = 0.25, acc_upper = 0.45)
# Usually automatic adjustment will work okay so long as burnin and n_draws both long enough
#Sample from posterior predictive distribution 
predict(run) <- predict(run, horizon = 20, conf_bands = c(0.05, 0.16))

## Can plot using following commands, but I will use ggplot for compatibility with forecast libraries
#plot(predict(run),area=TRUE,t_back = 190)

## With library(BVARverse) also have plotting option
#bv_ggplot(predict(run),t_back = 190)
#Alternative approach: BMR library
# Available at https://www.kthohr.com/bmr.html
#Uses slightly different formulation of and notation for Minnesota prior
#There are some useful options available in this library, 
# but it is not on CRAN and is not compatible with Kaggle's backend, so not run by default

#Convert to a data frame
bvarmacrodata <- data.matrix(macrodata) 

#Set up Minnesota-prior BVAR object, and sample from posterior by MCMC
# See https://www.kthohr.com/bmr_docs_vars_bvarm.html for syntax documentation: manual is out of date
bvar_obj <- new(bvarm)

#Construct BVAR with nlags lags and a constant
bvar_obj$build(data_endog=bvarmacrodata,cons_term=TRUE,p=nlags)

#Set random walk prior mean for all variables
coef_prior=c(1,1,1,1,1,1,1) 
# Set prior parameters (1,0.2,1) with harmonic decay
bvar_obj$prior(coef_prior=coef_prior,var_type=1,decay_type=1,HP_1=1,HP_2=0.2,HP_3=1,HP_4=2)

#Sample from BVAR with 10000 draws of Gibbs Sampler
bvar_obj$gibbs(10000)

#Construct BVAR Forecasts
bvarfcst<-BMR::forecast(bvar_obj,periods=20,shocks=TRUE,plot=FALSE,varnames=colnames(macrodata),percentiles=c(.05,.50,.95),
    use_mean=FALSE,back_data=0,save=FALSE,height=13,width=11)

#Warning: command is incredibly slow if plot=TRUE is on, fast otherwise
# Appears to be issue with plotting code in BMR package, which slows down plotting to visualize
# With too many lags and series, have to wait through hundreds of forced pauses

for (i in 1:nseries){
  BVAR<-ts(bvarfcst$forecast_mean[,i],start=c(2020,2),frequency=4,names=Series[i]) #Mean
  lcband<-ts(bvarfcst$plot_vals[,1,i],start=c(2020,2),frequency=4,names=Series[i]) #5% Lower confidence band
  ucband<-ts(bvarfcst$plot_vals[,3,i],start=c(2020,2),frequency=4,names=Series[i]) #95% Upper confidence band
  fdate<-time(lcband) #Extract date so geom_ribbon()knows what x value is
  bands<-data.frame(fdate,lcband,ucband) #Collect in data frame
}
#Construct Forecasts of Each Series by Univariate AR models, with 95% confidence intervals
ARIfcsts<-list()
for (i in 1:nseries) {
  ARIfcsts[[i]]<-forecast::forecast(ARIstatmodels[[i]],h=20,level=95)
}

# ADD VAR and ARI forecasts to plot

forecastseriesplots<-list()

#Plot ARI model forecast along with BVAR forecasts, plus respective 95% intervals
#Commented out alternative prior method
for (i in 1:nseries){
  BVAR<-ts(run$fcast$quants[3,,i],start=c(2020,2),frequency=4,names=Series[i]) #Mean
  lcband<-ts(run$fcast$quants[1,,i],start=c(2020,2),frequency=4,names=Series[i]) #5% Lower confidence band
  ucband<-ts(run$fcast$quants[5,,i],start=c(2020,2),frequency=4,names=Series[i]) #95% Upper confidence band
  fdate<-time(lcband) #Extract date so geom_ribbon() knows what x value is
  bands<-data.frame(fdate,lcband,ucband) #Collect in data frame
  #BMR_BVAR<-ts(bvarfcst$forecast_mean[,i],start=c(2020,2),frequency=4,names=Series[i]) #Mean
  #BMR_lcband<-ts(bvarfcst$plot_vals[,1,i],start=c(2020,2),frequency=4,names=Series[i]) #5% Lower confidence band
  #BMR_ucband<-ts(bvarfcst$plot_vals[,3,i],start=c(2020,2),frequency=4,names=Series[i]) #95% Upper confidence band
  #BMR_bands<-data.frame(fdate,BMR_lcband,BMR_ucband) #Collect in data frame
forecastseriesplots[[i]]<-autoplot(ARIfcsts[[i]])+autolayer(BVAR)+
  geom_ribbon(aes(x=fdate,ymin=lcband,ymax=ucband),data=bands,alpha=0.2,colour="red")+
  #autolayer(BMR_BVAR)+
  #geom_ribbon(aes(x=fdate,ymin=BMR_lcband,ymax=BMR_ucband),data=BMR_bands,alpha=0.2,colour="green")+
  labs(x="Date",y=colnames(macrodata)[i],title=Series[i])
}

grid.arrange(grobs=forecastseriesplots,nrow=4,ncol=2)

Interpretation

Conclusions

References