· 6 years ago · Jan 10, 2020, 08:00 PM
1### Recession indicator using unemployment. I take not-seasonally-adjusted
2### raw unemployment data (UNRATENSA) from FRED, use a MSTL decomposition to
3### decompose into Trend/Seasonality/Remainder, and look at the d(Trend)/dt,
4### the "velocity" of unemployment.
5###
6### For this indicator, the important feature *isn't* the peak,
7### but rather the point where unemployment velocity turns positive.
8### When the velocity passes through 0, that is a fairly reliable sign
9### that a recession is coming a few months down the road. In the
10### entire UNRATENSA dataset, there are only two "head fakes", which
11### can be mitigated by waiting for d(Trend)/dt to rise a small threshold
12### above zero.
13###
14### I have found that the "seasonally adjusted" data [UNRATE] contains
15### contains significant residual seasonality, so I use mstl() on the NSA data
16### [UNRATENSA] to ensure I have filtered out as much seasonality as possible.
17###
18### written by /u/MetricT (PM me for my actual email)
19
20### Set your FRED API key here. You may request an API key at:
21### https://research.stlouisfed.org/useraccount/apikeys
22api_key_fred <- "WHATEVER_YOUR_FRED_API_KEY_IS"
23
24####################################################################
25### Load necessary R packages, set the FRED API key, set start/end
26### dates for graph, add data frame with recession dates so we can
27### insert recession bars in our graph.
28####################################################################
29
30### We need the following packages for this example.
31packages <- c("fredr", "lubridate", "fredr", "ggplot2", "forecast",
32 "ggthemes", "tsibble", "dplyr", "magrittr", "broom", "scales")
33
34### Install packages if needed, then load them quietly
35new_packages <- packages[!(packages %in% installed.packages()[, "Package"])]
36if (length(new_packages)) install.packages(new_packages, quiet = TRUE)
37invisible(lapply(packages, "library",
38 quietly = TRUE,
39 character.only = TRUE,
40 warn.conflicts = FALSE))
41
42fredr_set_key(api_key_fred)
43
44### Specifiy the start/end date of the graph
45date_start <- as.Date("1948-01-01") # 1948-01-01 is earliest date of UNRATE
46date_end <- as.Date(now())
47
48####################################################################
49### Add recession bars to ggplot graphs
50####################################################################
51recessions_tibble <- tibble(
52
53 peak = as.Date(
54 c("1857-06-01", "1860-10-01", "1865-04-01", "1869-06-01",
55 "1873-10-01", "1882-03-01", "1887-03-01", "1890-07-01",
56 "1893-01-01", "1895-12-01", "1899-06-01", "1902-09-01",
57 "1907-05-01", "1910-01-01", "1913-01-01", "1918-08-01",
58 "1920-01-01", "1923-05-01", "1926-10-01", "1929-08-01",
59 "1937-05-01", "1945-02-01", "1948-11-01", "1953-07-01",
60 "1957-08-01", "1960-04-01", "1969-12-01", "1973-11-01",
61 "1980-01-01", "1981-07-01", "1990-07-01", "2001-03-01",
62 "2007-12-01")),
63
64 trough = as.Date(
65 c("1858-12-01", "1861-06-01", "1867-12-01", "1870-12-01",
66 "1879-03-01", "1885-05-01", "1888-04-01", "1891-05-01",
67 "1894-06-01", "1897-06-01", "1900-12-01", "1904-08-01",
68 "1908-06-01", "1912-01-01", "1914-12-01", "1919-03-01",
69 "1921-07-01", "1924-07-01", "1927-11-01", "1933-03-01",
70 "1938-06-01", "1945-10-01", "1949-10-01", "1954-05-01",
71 "1958-04-01", "1961-02-01", "1970-11-01", "1975-03-01",
72 "1980-07-01", "1982-11-01", "1991-03-01", "2001-11-01",
73 "2009-06-01")
74 )
75 )
76
77geom_recession_bars <- function (date_df_start, date_df_end) {
78
79 date_df_start <- as.Date(date_df_start, origin="1970-01-01")
80 date_df_end <- as.Date(date_df_end, origin="1970-01-01")
81
82 recessions_df <- read.table(textConnection(
83 "peak, trough
84 1857-06-01, 1858-12-01
85 1860-10-01, 1861-06-01
86 1865-04-01, 1867-12-01
87 1869-06-01, 1870-12-01
88 1873-10-01, 1879-03-01
89 1882-03-01, 1885-05-01
90 1887-03-01, 1888-04-01
91 1890-07-01, 1891-05-01
92 1893-01-01, 1894-06-01
93 1895-12-01, 1897-06-01
94 1899-06-01, 1900-12-01
95 1902-09-01, 1904-08-01
96 1907-05-01, 1908-06-01
97 1910-01-01, 1912-01-01
98 1913-01-01, 1914-12-01
99 1918-08-01, 1919-03-01
100 1920-01-01, 1921-07-01
101 1923-05-01, 1924-07-01
102 1926-10-01, 1927-11-01
103 1929-08-01, 1933-03-01
104 1937-05-01, 1938-06-01
105 1945-02-01, 1945-10-01
106 1948-11-01, 1949-10-01
107 1953-07-01, 1954-05-01
108 1957-08-01, 1958-04-01
109 1960-04-01, 1961-02-01
110 1969-12-01, 1970-11-01
111 1973-11-01, 1975-03-01
112 1980-01-01, 1980-07-01
113 1981-07-01, 1982-11-01
114 1990-07-01, 1991-03-01
115 2001-03-01, 2001-11-01
116 2007-12-01, 2009-06-01"),
117 sep = ",",
118 colClasses = c("Date", "Date"),
119 header = TRUE)
120
121 recessions_trim <- subset(recessions_df,
122 peak >= min(date_df_start) &
123 trough <= max(date_df_end))
124
125 if (nrow(recessions_trim) > 0) {
126 recession_bars = geom_rect(data = recessions_trim, inherit.aes = F,
127 fill = "darkgray", alpha = 0.25,
128 aes(xmin = peak, xmax = trough, ymin = -Inf, ymax = +Inf))
129 } else {
130 recession_bars = geom_blank()
131 }
132}
133
134####################################################################
135### Fetch unemployment data from FRED
136####################################################################
137
138### Load unemployment from FRED
139target <- c("UNRATENSA", "US")
140
141series <- target[1]
142name <- target[2]
143
144data <- fredr(series_id = series,
145 frequency = "m",
146 observation_start = date_start,
147 observation_end = date_end)
148
149date <- data %>% as_tsibble(index = "date") %>% pull("date")
150values <- data %>% as_tsibble(index = "date") %>% pull("value")
151
152# Decompose the data into trend/seasonal/remainder
153values_decomposed <- ts(values, frequency = 12,
154 start = c(year(min(date)), month(min(date)))) %>%
155 mstl()
156
157# Get the seasonal, trend, and remainder components
158seasonal <- values_decomposed %>% seasonal()
159trend <- values_decomposed %>% trendcycle()
160remainder <- values_decomposed %>% remainder()
161
162unemployment <- trend
163
164####################################################################
165### Calculations to derive unemployment "velocity"
166####################################################################
167diff_trend <- unemployment %>% diff()
168diff_date <- date %>% tail(n = length(diff_trend))
169
170### Unemployment is reported with units of 0.1%. ie, unemployment is 3.5%, 3.4%,
171### etc. As a consequence, that quantization introduces "stairstepping" when
172### taking the derivative. As a first pass, I use Friedman's SuperSmoother
173### to eliminate it. There's almost certainly a more mathematically correct
174### way to do this, but damn it Jim, I'm a physicist, not an economist...
175### Just comment the two lines below out if you want use raw data instead of the
176### smoothed data. It works the same and arrives at the same conclusions, but
177### the unfiltered data yields slightly uglier graphs.
178#tmp <- supsmu(decimal_date(diff_date), diff_trend, span = 0.0085)
179tmp <- supsmu(decimal_date(diff_date), diff_trend, span = 3/length(diff_trend))
180diff_trend <- tmp$y
181
182velocity_df <- data.frame(diff_date, diff_trend)
183
184####################################################################
185### Find locations where the unemployment velocity line grows
186### through x = 0
187####################################################################
188axis_crossing_locations <- function(date, value) {
189
190 # Set a threshold to filter out minor incursions about x = 0
191 # This can improve the accuracy of the indicator at the cost
192 # of reducing the time between the predictor triggering and the
193 # recession date. As my original reason for writing this was
194 # trying to make sure I get my money out of the market before
195 # the stock market starts sliding downhill, I personally use
196 # a threshold of 0 and take a chance that it could be wrong.
197 # Anyone more concerned with accuracy should use a higher
198 # threshold
199 # threshold <- 0
200
201 # The threshold necessary to eliminate the two false positives from the 1960's
202 threshold <- 0.017
203
204 # Shift the value by the threshold
205 adj_value <- value - threshold
206
207 # See if the sign changes before/after, which indicates it has
208 # passed through x = 0
209 sign <- sign(value)
210 adj_sign <- sign(adj_value)
211
212 # Computer time derivative of sign changes
213 diff_sign <- diff(sign)
214 diff_adj_sign <- diff(adj_sign)
215
216 # This will pick up *growing* unemployment. To pick up
217 # *falling* unemployment, use -2. To find both, use
218 # abs(diff_sign). We want to add one month to the
219 # date this gives us, thus the "%m+% months(1)" bit.
220 locs <- date[diff_adj_sign == 2] %m+% months(1)
221 locs
222}
223
224axiscrossings <- as.Date(axis_crossing_locations(diff_date, diff_trend))
225
226
227####################################################################
228### Graph: Trend Unemployment Velocity
229####################################################################
230c1 <- paste("U.S. Bureau of Labor Statistics, Unemployment Rate in ",
231 name, " [", series, "]\n", sep = "")
232c2 <- "Retrieved from FRED, Federal Reserve Bank of St. Louis;\n"
233c3 <- paste("https://fred.stlouisfed.org/series/", series, "\n")
234c4 <- paste("Data retrieved on", format(Sys.time(), "%B %d, %Y at %I:%M %p %Z"))
235
236s1 <- "Indicator (vertical dotted line) shows imminent recession when velocity rises above 0\n"
237s2 <- "The farther above/below 0, the faster unemployment is growing/falling\n"
238s3 <- "Recessions marked with vertical bars\n"
239
240title <- paste(name, "Trend Unemployment Rate Velocity")
241#title <- "Step 3: Instead of crossing axis at x = 0, add a small threshold\n (~0.032) to eliminate false positives in 1960's and 1995"
242subtitle <- paste(s1, s2, s3, sep = "")
243#subtitle <- s3
244xlab <- "Year"
245ylab <- "Percent/Month"
246caption <- paste(c1, c2, c3, c4, sep = "")
247
248p <- ggplot(velocity_df, aes(x = diff_date, y = diff_trend)) +
249
250 ### Plot unemployment velocity
251 geom_line(size = 1.3, color = "red3") +
252
253
254 ### Add a line at velocity = 0. Above this line, unemployment is rising, and
255 ### the farther above the line, the faster it rises. Similarly, below the
256 ### line unemployment is falling, and the further below the line, the faster
257 ### it falls
258 geom_hline(yintercept = 0, size = 1) +
259
260 ### Draw the recession indicators when the unemployment velocity has
261 ### has grown beyond x = 0
262 geom_vline(xintercept = axiscrossings, linetype = "dashed") +
263
264 # Add recession bars
265 geom_recession_bars(min(velocity_df$diff_date), max(velocity_df$diff_date)) +
266
267 ### Misc graph stuff
268 theme_economist() +
269 scale_x_date(breaks = pretty_breaks(10), limits = c(date_start, date_end)) +
270 labs(title = title, subtitle = subtitle, caption = caption,
271 x = xlab, y = ylab)
272
273print(p)