Automatically plotting and forecasting any BLS series
William Campbell
2024-07-26
This document showcases a pipeline to automatically forecast and plot any data series available from the Bureau of Labor Statistics (BLS). The pipeline is built on top of BLS API’s, as well as the forecast and ggplot2 packages in R. Using a BLS API key is highly recommended to avoid query limitations.
arima_plots function
The arima_plots function allows for generation of an auto-fitted ARIMA model which is plotted along with historical data for any BLS series. Various arguments are accepted to customize the returned plot. Source code for the first version is available in the Appendix.
In future versions support for FRED data, smoothing models, and other chart types will be added. Additionally, support for simple and advanced visualization without forecasting is under development.
#Example usage
arima_plots(series_names = c('LNS14000000','CUUR0000SA0','JTS000000000000000JOL','WPUFD49207'),
titles =c('Unemployment Rate (Seasonally Adjusted)','CPI for All Urban Consumers',
'JOLTS (Seasonally Adjusted)','PPI Finished Goods'),
include = 40,
mom_growth_rate = c(FALSE,TRUE,FALSE,TRUE),
yoy_growth_rate = c(FALSE,FALSE,FALSE,FALSE),
ylabel = c('Percent','Percent Change MoM','Number in Thousands','Percent Change MoM'))
The arima_plots() function also allows year over year and monthly growth rates to be plotted as shown below:
#Plotting CPI three different ways
arima_plots(series_names = c('CUUR0000SA0','CUUR0000SA0','CUUR0000SA0'),
titles = c('CPI Index 1982-84 = 100','CPI Month over Month',
'CPI Year over Year'), ylabel = c('Index Value', 'Growth Rate', 'Growth Rate'),
include = 40,
mom_growth_rate = c(FALSE,TRUE,FALSE),
yoy_growth_rate = c(FALSE,FALSE,TRUE),
n_col = 1,
summary = TRUE)
| Series | CurrentValue | OneYearAvg | OneYearHigh | OneYearLow | Forecast1 | Forecast2 | Forecast3 |
|---|---|---|---|---|---|---|---|
| CPI Index 1982-84 = 100 | 314.18 | 309.57 | 314.18 | 305.69 | 313.94 | 314.34 | 314.6 |
| CPI Month over Month | 0.03 | 0.24 | 0.64 | -0.2 | -0.09 | 0.11 | 0.07 |
| CPI Year over Year | 2.97 | 3.3 | 3.7 | 2.97 | 2.63 | 2.14 | 1.86 |
The arima_plots function makes creating economic outlook reports simple. The following example shows how one can generate such reports with only three lines of code.
Economic Outlook Dashboard, July 2024
Employment
| Series | CurrentValue | OneYearAvg | OneYearHigh | OneYearLow | Forecast1 | Forecast2 | Forecast3 |
|---|---|---|---|---|---|---|---|
| Unemployment Rate (Seasonally Adjusted) | 4.1 | 3.81 | 4.1 | 3.5 | 4.21 | 4.3 | 4.38 |
| JOLTS (Seasonally Adjusted) | 8140 | 8756.25 | 9358 | 7919 | 8139.09 | 8039.92 | 8093.51 |
Productivity
| Series | CurrentValue | OneYearAvg | OneYearHigh | OneYearLow | Forecast1 | Forecast2 | Forecast3 |
|---|---|---|---|---|---|---|---|
| Output Per Hour - Non-farm Business Productivity | 0.2 | 0.33 | 4.6 | -5.7 | 1.59 | 1.59 | 1.59 |
| Nonfarm Business Unit Labor Costs | 4 | 3.92 | 8.3 | -2.8 | 1.14 | 2.02 | 1.75 |
Inflation
| Series | CurrentValue | OneYearAvg | OneYearHigh | OneYearLow | Forecast1 | Forecast2 | Forecast3 |
|---|---|---|---|---|---|---|---|
| CPI for All Urban Consumers | 2.97 | 3.3 | 3.7 | 2.97 | 2.63 | 2.14 | 1.86 |
| PPI Final Demand (Seasonally Adjusted) | 2.65 | 1.63 | 2.65 | 0.77 | 2.18 | 1.71 | 1.59 |
Appendix
arima_plots = function(series_names,
titles = '',
autofit = TRUE, ### not yet available
seasonal = FALSE, ### not yet available
n_ahead = 8,
include = 20,
frequency = 'month',
ylabel = '',
mom_growth_rate=c(FALSE),
yoy_growth_rate = c(FALSE),
years_included = 9,
n_col=2,
summary = FALSE
){
library(blsR)
library(forecast)
library(ggplot2)
library(zoo)
library(gridExtra)
library(kableExtra)
plot_list <<- list()
summary_stats = data.frame(
Series = character(length(series_names)),
CurrentValue = numeric(length(series_names)),
OneYearAvg = numeric(length(series_names)),
OneYearHigh = numeric(length(series_names)),
OneYearLow = numeric(length(series_names)),
Forecast1 = numeric(length(series_names)),
Forecast2 = numeric(length(series_names)),
Forecast3 = numeric(length(series_names)),
stringsAsFactors = FALSE
)
for (i in seq_along(series_names)){
#manage start and end periods
current_year = as.numeric(format(Sys.Date(), "%Y"))
start_year = current_year - years_included
#query series from BLS
#querying one series at a time may not be the most efficient, however users without
#an API key are limited to one series query at a time
query = query_series(series_names[i],start_year = start_year, end_year = current_year)
data = bls_request(query)
data_as_table = data_as_tidy_table(data$series[[1]]$data)
#manage frequency for time series conversion
if (frequency == 'month') {
frequency_num = 12
} else if (frequency == 'quarter') {
frequency_num = 4
}
#convert to time series
series = ts(data_as_table$value, start = start_year,frequency = frequency_num)
#calculate growth rate if specified in argument
if (mom_growth_rate[i] ==TRUE){
series = diff(log(series))*100
series = na.omit(series)
}
if (yoy_growth_rate[i]==TRUE){
series = (series / lag(series, -12) - 1) * 100
series = na.omit(series)
}
#fit and forecast arima
fit = auto.arima(series)
forecast_data = forecast(fit, h=n_ahead)
#prepare data for plotting - truncate to desired plotting length
subset_data = tail(series,include)
actual_data <- data.frame(
Date = as.Date(as.yearmon(time(subset_data))),
Value = as.numeric(subset_data)
)
#prepare dataframes for use with ggplot
last_date = as.Date(as.yearmon(time(subset_data)[length(subset_data)]))
forecast_dates = seq(last_date + 1, by = frequency, length.out = length(forecast_data$mean))
forecast_data_frame <- data.frame(
Date = forecast_dates,
Forecast = as.numeric(forecast_data$mean),
Lower80 = as.numeric(forecast_data$lower[,1]),
Upper80 = as.numeric(forecast_data$upper[,1]),
Lower95 = as.numeric(forecast_data$lower[,2]),
Upper95 = as.numeric(forecast_data$upper[,2])
)
plot =ggplot() +
geom_line(data = actual_data, aes(x = Date, y = Value), color = "black") +
geom_line(data = forecast_data_frame, aes(x = Date, y = Forecast), color = "red") +
geom_ribbon(data = forecast_data_frame, aes(x = Date, ymin = Lower80, ymax = Upper80), fill = "orange", alpha = 0.3) +
geom_ribbon(data = forecast_data_frame, aes(x = Date, ymin = Lower95, ymax = Upper95), fill = "yellow", alpha = 0.3) +
ggtitle(titles[i]) +
xlab("Year") +
ylab(ylabel[i]) +
theme_minimal()+
theme(
plot.background = element_rect(fill = "white", color = NA),
panel.background = element_rect(fill = "white", color = NA),
panel.grid.major = element_line(color = "gray95"),
panel.grid.minor = element_line(color = "gray95"),
plot.title = element_text(size = 16,hjust = .5),
plot.title.position = 'plot',
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
axis.text = element_text(size = 12)
)
plot_list[[i]] <<- plot ### send plots to a global list so user can reproduce individually
# Calculate summary statistics
one_year_data = tail(series, 12)
current_value = tail(series, 1)
one_year_avg = mean(one_year_data, na.rm = TRUE)
one_year_high = max(one_year_data, na.rm = TRUE)
one_year_low = min(one_year_data, na.rm = TRUE)
# Forecast values
forecast_1 = ifelse(length(forecast_data$mean) >= 1, forecast_data$mean[1], NA)
forecast_2 = ifelse(length(forecast_data$mean) >= 2, forecast_data$mean[2], NA)
forecast_3 = ifelse(length(forecast_data$mean) >= 3, forecast_data$mean[3], NA)
# Store summary statistics in the data frame
summary_stats[i, ] = c(titles[i],
round(current_value, 2),
round(one_year_avg, 2),
round(one_year_high, 2),
round(one_year_low, 2),
round(forecast_1, 2),
round(forecast_2, 2),
round(forecast_3, 2))
}
# Display the summary statistics table
library(knitr)
summary_stats[] <- lapply(summary_stats, as.character) # Convert all columns to character for proper display
#plot all plots
do.call(grid.arrange, c(plot_list, ncol = n_col))
if (summary==TRUE){
print(
kable(summary_stats, caption = "Summary Statistics for Each Series", align = 'c') %>%
kable_styling(full_width = FALSE, position = "center", bootstrap_options = c("striped", "hover", "condensed", "responsive"))
)
}
}