Loading packages, defining colors and using data

We first clear the workspace using rm(list = ls())and then include all packages we need. If a package is missing in your R distribution (which is quite likely initially), just use install.packages("package_name") with the respective package name to install it on your system. If you execute the code in the file install_packages.R, then all necessary packages will be installed into your R distribution. If the variable export_graphs is set to TRUE, then the graphs will be exported as pdf-files. In addition, we define a set of colors here to make graphs look more beautiful. Furthermore, we define some axis breaks for graphs to look nicer.

rm(list = ls())
library(zoo)
library(xts)
library(dynlm)
library(reshape2)
library(base)
library(ggplot2)
library(grid)
library(scales)
library(stringr)
library(tidyverse)
library(fredr)
library(mFilter)
library(formattable)
library(nortsTest)

# should graphs be exported to pdf
export_pdf <- FALSE

# define some colors
mygreen <- "#00BA38"
myblue  <- "#619CFF"
myred   <- "#F8766D"
  
# define plotting breaks
xbreaks <- c(seq(from = as.Date("1950-01-01"), to = as.Date("2020-01-01"),by = "10 years"))

Pull data from FRED database

In a first step, we need to pull relevant data from the FRED database. We again do this using the function get_fred_data, which we already discussed in Chapter 6. This function allows us to pull multiple data series at the same time as long as they have the same starting and end date and the same frequency. We first pull a series of quarterly data consisting of nominal GDP and its components, the deflators for different GDP components and some employment and labor compensation statistics.

source("chap06_func_get_fred_data.R")

########### 
# Quarterly time series
###########

# nominal time series
selection <- data.frame(series = "GDP", names = "GDP")
selection <- rbind(selection, c("PCEC", "CONS"))
selection <- rbind(selection, c("PCND", "CONS_N"))
selection <- rbind(selection, c("PCESV", "CONS_S"))
selection <- rbind(selection, c("GPDI", "INV"))
selection <- rbind(selection, c("FPI", "INV_F"))
selection <- rbind(selection, c("GCE", "GOV"))
selection <- rbind(selection, c("TB3MS", "NOMRATE"))
selection <- rbind(selection, c("EXPGS", "EXP"))
selection <- rbind(selection, c("IMPGS", "IMP"))

# price deflators
selection <- rbind(selection, c("GDPDEF", "P_GDP"))
selection <- rbind(selection, c("DPCERD3Q086SBEA", "P_CONS"))
selection <- rbind(selection, c("DNDGRD3Q086SBEA", "P_CONS_N"))
selection <- rbind(selection, c("DSERRD3Q086SBEA", "P_CONS_S"))
selection <- rbind(selection, c("A006RD3Q086SBEA", "P_INV"))
selection <- rbind(selection, c("A007RD3Q086SBEA", "P_INV_F"))
selection <- rbind(selection, c("A822RD3Q086SBEA", "P_GOV"))
selection <- rbind(selection, c("A020RD3Q086SBEA", "P_EXP"))
selection <- rbind(selection, c("A021RD3Q086SBEA", "P_IMP"))

# employment and labor compensation
selection <- rbind(selection, c("HOABS", "HOURS"))
selection <- rbind(selection, c("HOANBS", "HOURS_NF"))
selection <- rbind(selection, c("A358RX1Q020SBEA", "VA_NF"))
selection <- rbind(selection, c("A439RC1Q027SBEA", "VA_NET"))
selection <- rbind(selection, c("A442RC1Q027SBEA", "EMPCOMP"))

# pull all quarterly data series from FRED
fred_q <- get_fred_data(selection$series, selection$names, "1948-01-01", "2023-12-31", "q")
## [1] "Loading series PCEC"
## [1] "Loading series PCND"
## [1] "Loading series PCESV"
## [1] "Loading series GPDI"
## [1] "Loading series FPI"
## [1] "Loading series GCE"
## [1] "Loading series TB3MS"
## [1] "Loading series EXPGS"
## [1] "Loading series IMPGS"
## [1] "Loading series GDPDEF"
## [1] "Loading series DPCERD3Q086SBEA"
## [1] "Loading series DNDGRD3Q086SBEA"
## [1] "Loading series DSERRD3Q086SBEA"
## [1] "Loading series A006RD3Q086SBEA"
## [1] "Loading series A007RD3Q086SBEA"
## [1] "Loading series A822RD3Q086SBEA"
## [1] "Loading series A020RD3Q086SBEA"
## [1] "Loading series A021RD3Q086SBEA"
## [1] "Loading series HOABS"
## [1] "Loading series HOANBS"
## [1] "Loading series A358RX1Q020SBEA"
## [1] "Loading series A439RC1Q027SBEA"
## [1] "Loading series A442RC1Q027SBEA"


Since data on the capital stock is only available at the annual frequency, we pull some data on aggregate production also at the annual frequency level.

########### 
# Annual time series
###########

selection <- data.frame(series = "GDP", names = "GDP")
selection <- rbind(selection, c("K1TTOTL1ES000", "CAPSTOCK"))
selection <- rbind(selection, c("B4701C0A222NBEA", "TOTALHOURS"))
selection <- rbind(selection, c("PAYEMS", "TOTALEMPS"))
selection <- rbind(selection, c("A442RC1A027NBEA", "COMPENSATION"))

# pull all annual data series from FRED
fred_a <- get_fred_data(selection$series, selection$names, "1950-01-01", "2022-12-31", "a")
## [1] "Loading series K1TTOTL1ES000"
## [1] "Loading series B4701C0A222NBEA"
## [1] "Loading series PAYEMS"
## [1] "Loading series A442RC1A027NBEA"
# calculate hours per employee
fred_a$HOURS = fred_a$TOTALHOURS/fred_a$TOTALEMPS*1000

Deflating variables to real values

We first deflate all components of GDP to get to their real values. For them to sum up to total real GDP, we again need to make a slight adjustment. We already discussed this in Chapter 6.

# calculate real variables by deflating nominal ones
for(var in c("GDP", "CONS", "CONS_N", "CONS_S", "INV", "INV_F", "GOV", "EXP", "IMP")) {
  fred_q[paste("R_", var, sep="")] = fred_q[var]/fred_q[paste("P_", var, sep="")]*100
}

# calculate adjusted real GDP (because of different deflators that don't finally add up correctly)
factor <- fred_q$R_GDP/(fred_q$R_CONS +fred_q$R_INV + fred_q$R_GOV + fred_q$R_EXP - fred_q$R_IMP)

# scale series
for(var in c("GDP", "CONS", "INV", "GOV", "EXP", "IMP")) {
  fred_q[paste("R_", var, sep="")] = fred_q[paste("R_", var, sep="")]*factor
}

The cyclical properties of government consumption

We now study the cyclical properties of (the log of) government consumption. To this end, we first extract the cycle from the time series of government spending using the HP filter. We then plot the cyclical component of government spending over time. Government spending appears to have some persistence. However, there seems to be some excess motion in government spending in the first 40 quarters or so (right after World War II). For the following analysis, we hence only use the time series for government spending starting at quarter 40.

# filter government spending
ts <- xts(x = log(fred_q$R_GOV), order.by = as.Date(fred_q$date))
hpf <- hpfilter(ts, freq=1600, type="lambda")
fred_q$gov_cycle <- hpf$cycle
gov_exp <- hpf$cycle[40:length(hpf$cycle)]

# Plot cycle of government expenditure
myplot <- ggplot(data = fred_q) + 
  geom_vline(xintercept=40, color="gray", linewidth=0.5) +
  geom_hline(yintercept=0, color="gray", linewidth=0.5) +
  geom_line(aes(x=quarter, y=gov_cycle), color="darkblue", linewidth=1) +
  labs(x = "Quarter t",
       y = "Cycle of log(G_t)") +
  coord_cartesian(ylim=c(-0.15, 0.15)) +
  scale_y_continuous(breaks=seq(-0.15, 0.15, 0.05), label=comma) +  
  theme_bw()

# print the plot
print(myplot)


In order to trace out the persistence of the cyclical component of government spending, we calculate its partial autocorrelations (for the data starting at quarter 40). As we can see, the partial autocorrelation is large and significant at a lag of one period. All other autocorrelations are minor and (if at all) hardly significant. Hence, assuming an AR(1) process for government spending seems to be safe.

# generate partial autocorrelation plot
myplot <- ggpacf(gov_exp[40:length(gov_exp)]) +
  coord_cartesian(ylim=c(-0.5, 1)) +  
  scale_y_continuous(breaks=seq(-0.5, 1, 0.25)) +  
  labs(x = "Lag",
       y = "Partial Auto-Correlation") +
  theme_bw()
print(myplot)

Calculate calibration targets

Now we want to use our FRED data to calculate calibration targets. For most of it, we just use the long-run averages (mean). In addition, we can calculate the persistence parameter of the cyclical components of the log of government consumption with a simple regression using the lagged value (lag) to explain the current value. We finally show all the calibration targets in a table.

# calculate relevant statistics
dat <- array(dim=c(8))

# consumption share in GDP
dat[1] <- mean(fred_q$CONS/fred_q$GDP)
dat[2] <- mean(fred_q$INV/fred_q$GDP)
dat[3] <- mean(fred_q$GOV[40:length(fred_q$GOV)]/fred_q$GDP[40:length(fred_q$GDP)])
dat[4] <- mean(fred_a$HOURS)/(365*(24-8))
dat[5] <- mean(fred_q$EMPCOMP/fred_q$VA_NET)

# attention: capital stock measured in millions and GDP in billions
dat[6] <- mean(fred_a$CAPSTOCK/(fred_a$GDP/4))/1000

# calculate autocorrelation of government spending
reg <- lm(gov_exp[40:length(gov_exp)] ~ lag(gov_exp[40:length(gov_exp)]))
dat[7] <- reg$coefficients[2]
dat[8] <- sd(gov_exp[40:length(gov_exp)])*sqrt(1 - dat[7]^2)


# generate data table
table_data <- data.frame("Statistics" = 
                           c("Consumption C/Y", "Investment I/Y", 
                             "Government purchases G/Y", "Hours worked (share of endowment)",
                             "Labor share corporate sector wL/Y", "Capital K/Y (quarterly)",
                             "Government purchases: autocorrelation rho_G", "Government purchases: sd sigma_G"))
table_data <- cbind(table_data, "US Data" = digits(dat, 4, format="f"))


# output table
formattable(table_data, align= c("l", "r"))
Statistics US Data
Consumption C/Y 0.6357
Investment I/Y 0.1739
Government purchases G/Y 0.2037
Hours worked (share of endowment) 0.3222
Labor share corporate sector wL/Y 0.7030
Capital K/Y (quarterly) 11.6760
Government purchases: autocorrelation rho_G 0.8415
Government purchases: sd sigma_G 0.0073

Calculate calibrated parameters

From the data calculated above, we need to derive the model parameters. All parameters are directly linked to one or more moments. We discuss these linkages in the lecture. Having calculated all parameters, we print them to a table.

# share or labor compensation -> production share
alpha <- 1 - dat[5]

# investment share in GDP -> depreciation
delta <- dat[2]/dat[6]

# capital share in production -> time preference rate
theta <- alpha/dat[6] - delta

# time discount factor
beta  <- 1/(1+theta) 

# labor hours share -> disutility of labor
chi   <- (1-alpha)/(1 - dat[2] - dat[3]) *(1-dat[4])/dat[4]

# government share in production
g_y   <- dat[3]

# government expenditure/GDP -> G_bar
G_bar <- log(dat[3]*dat[6]^(alpha/(1-alpha))*dat[4])

# autocorrelation parameter government expenditure
rho_G    <- dat[7]

# variance of government expenditure
sigma_G  <- dat[8]

# autocorrelation parameter of technology shocks
rho_A    <- 0.95


# generate calibrated value table
table_data <- data.frame("Parameter" = 
                           c("Production parameter (alpha)", "Depreciation rate (delta)", 
                             "Time discount factor (beta)", "Disutility of labor (chi)",
                             "Government share (g_y)", "Mean gov. expend. process (G_bar)",
                             "Autocorrelation gov. expend. rho_G", "Standard deviation gov. expend. sigma_G"))
table_data <- cbind(table_data, "Value" = digits(c(alpha, delta, beta, chi, g_y, G_bar, rho_G, sigma_G), 4, format="f"))


# output table
formattable(table_data, align= c("l", "r"))
Parameter Value
Production parameter (alpha) 0.2970
Depreciation rate (delta) 0.0149
Time discount factor (beta) 0.9896
Disutility of labor (chi) 2.3763
Government share (g_y) 0.2037
Mean gov. expend. process (G_bar) -1.6858
Autocorrelation gov. expend. rho_G 0.8415
Standard deviation gov. expend. sigma_G 0.0073

Calculate the steady state

Using the calibrated parameters, we can easily derive the non-stochastic steady state of the RBC model. This steady state has a closed form solution.

# interest rate
r  <- theta

# capital intensity
KY <- alpha/(theta + delta)

# investment share
IY <- delta*KY

# government share
GY <- dat[3]

# consumption share
CY <- 1 - IY - GY

# labor supply
L  <- (1-alpha)/(1-alpha + chi*CY)

# output level
Y  <- KY^(alpha/(1-alpha)) * L

# wage rate
w = (1-alpha)*Y/L

# share of labor compensation
wL <- w*L/Y


# generate data table
table_data <- data.frame("Statistics" = 
                           c("Consumption C/Y", "Investment I/Y", 
                             "Government purchases G/Y", "Hours worked (share of endowment)",
                             "Labor share corporate sector wL/Y", "Capital K/Y (quarterly)",
                             "Interest rate (quarterly)"))

table_data <- cbind(table_data, "Model Value" = digits(c(CY, IY, GY, L, wL, KY, r), 4, format="f"))

# output table
formattable(table_data, align= c("l", "r"))
Statistics Model Value
Consumption C/Y 0.6225
Investment I/Y 0.1739
Government purchases G/Y 0.2037
Hours worked (share of endowment) 0.3222
Labor share corporate sector wL/Y 0.7030
Capital K/Y (quarterly) 11.6760
Interest rate (quarterly) 0.0105

Save parameters and steady state values to file

To work with our parameter set as well as the calculated steady state values throughout this entire chapter, we save the variables to an .RData file. We can use the function save to do this. In this function, we can specify a filename as well as a list of variable values that should be saved in this file.

save(file="chap07_parameters.RData", 
     list=c("alpha", "delta", "beta", "chi", "g_y", "G_bar", "rho_G", "sigma_G", "rho_A",
            "CY", "IY", "GY", "L", "r", "w", "KY", "Y"))

Filter time series using HP filter to replicate Hansen and Wright (1992, Tab. 2)

Last but not least, we compute some additional business cycle statistics which later on will allow us to compare our real business cycle model with the data. To this end, we first use an HP filter to filter the time series of GDP, consumption, investment, employment, and productivity.

# filter GDP
ts <- xts(x = log(fred_q$R_GDP), order.by = as.Date(fred_q$date))
hpf <- hpfilter(ts, freq=1600, type="lambda")
fred_q$GDP_cycle <- hpf$cycle 

# filter consumption
ts <- xts(x = log(fred_q$R_CONS_N + fred_q$R_CONS_S), order.by = as.Date(fred_q$date))
hpf <- hpfilter(ts, freq=1600, type="lambda")
fred_q$CONS_cycle <- hpf$cycle 

# filter investment
ts <- xts(x = log(fred_q$R_INV_F), order.by = as.Date(fred_q$date))
hpf <- hpfilter(ts, freq=1600, type="lambda")
fred_q$INV_cycle <- hpf$cycle 

# filter employment in non farm business
ts <- xts(x = log(fred_q$HOURS_NF), order.by = as.Date(fred_q$date))
hpf <- hpfilter(ts, freq=1600, type="lambda")
fred_q$HOURS_cycle <- hpf$cycle 

# filter productivity in non-farm business
ts <- xts(x = log(fred_q$VA_NF/fred_q$HOURS_NF), order.by = as.Date(fred_q$date))
hpf <- hpfilter(ts, freq=1600, type="lambda")
fred_q$PROD_cycle <- hpf$cycle


With the filtered components, we can calculate standard deviations of all variables relative to the standard deviation of GDP. This allows us to compare the volatility of different time series generated by the data and the model. In addition, we calculate the correlation of all variables with GDP. We finally print our results to a table.

# calculate relevant statistics
emp_data <- array(dim=c(6, 2))

# consumption share in GDP
emp_data[1, 1] <- sd(fred_q$GDP_cycle)
emp_data[2, 1] <- sd(fred_q$CONS_cycle) /sd(fred_q$GDP_cycle)
emp_data[3, 1] <- sd(fred_q$INV_cycle)  /sd(fred_q$GDP_cycle)
emp_data[4, 1] <- sd(fred_q$PROD_cycle) /sd(fred_q$GDP_cycle)
emp_data[5, 1] <- sd(fred_q$HOURS_cycle)/sd(fred_q$GDP_cycle)
emp_data[6, 1] <- sd(fred_q$HOURS_cycle) /sd(fred_q$PROD_cycle)

emp_data[1, 2] <- 1
emp_data[2, 2] <- cor(fred_q$CONS_cycle, fred_q$GDP_cycle)
emp_data[3, 2] <- cor(fred_q$INV_cycle , fred_q$GDP_cycle)
emp_data[4, 2] <- cor(fred_q$PROD_cycle, fred_q$GDP_cycle)
emp_data[5, 2] <- cor(fred_q$HOURS_cycle,fred_q$GDP_cycle)
emp_data[6, 2] <- cor(fred_q$HOURS_cycle,fred_q$PROD_cycle)


# generate data table
table_data <- data.frame("Statistics" = 
                           c("sd(Y) in %", 
                             "sd(C)/sd(Y)", "corr(C, Y)",
                             "sd(I)/sd(Y)", "corr(I, Y)",
                             "sd(w)/sd(Y)", "corr(w, Y)",
                             "sd(L)/sd(Y)", "corr(L, Y)",
                             "sd(L)/sd(w)", "corr(L, w)"))

table_data <- cbind(table_data, "US Data" = digits(c(emp_data[1, 1]*100, emp_data[2, 1],
                                emp_data[2, 2], emp_data[3, 1], emp_data[3, 2], emp_data[4, 1], emp_data[4, 2],
                                emp_data[5, 1], emp_data[5, 2], emp_data[6, 1], emp_data[6, 2]), 2, format="f"))


# output table
formattable(table_data, align= c("l", "r"))
Statistics US Data
sd(Y) in % 1.89
sd(C)/sd(Y) 0.56
corr(C, Y) 0.71
sd(I)/sd(Y) 2.47
corr(I, Y) 0.86
sd(w)/sd(Y) 0.60
corr(w, Y) 0.29
sd(L)/sd(Y) 1.11
corr(L, Y) 0.84
sd(L)/sd(w) 1.86
corr(L, w) -0.22


Finally, we save our results in a file named chap07_empirics.RData for later use.

# save empirical data
save(file="chap07_empirics.RData", 
     list=c("emp_data"))