Below, I discuss some `R`

code for illustrating how to implement the business cycle indicator proposed by de Carvalho, Rodrigues, and Rua (2012); for the interpretation of the quantities computed in this script, the reader is referred to the article. Before running the code chunks below, start by cleaning workspace and install the following packages (if not installed).

```
rm(list = ls())
if (!require("mFilter")) install.packages("mFilter")
```

The `mFilter`

package is here required to compare the CRR (Carvalho–Rodrigues–Rua) indicator with that of Christiano-Fitzgerald. For reproducibility reasons, I list below the information about `R`

, the OS, and loaded packages:

```
set.seed(1)
sessionInfo()
```

```
## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Sierra 10.12.1
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] mFilter_0.1-3
##
## loaded via a namespace (and not attached):
## [1] magrittr_1.5 assertthat_0.1 tools_3.3.2 htmltools_0.3.5
## [5] yaml_2.1.13 tibble_1.2 Rcpp_0.12.5 stringi_1.1.1
## [9] rmarkdown_1.1 knitr_1.15 stringr_1.0.0 digest_0.6.9
## [13] evaluate_0.10
```

In the code chunks below, I follow the 80 characters per line standard.

This is the key function to be used below:

```
ssa <- function(y, l, m) {
# SSA needs as input a time-series y
# SSA also needs two parameters to be defined by the user:
# - windows length l
# - number of leading eigentriples, m.
n <- length(y)
k <- n - l + 1
# Take I (Decomposition)
# Step 1 (embedding)
# constructs the trajectory matrix X, i.e. a lagged version of the
# original time-series y; this is a Hankel matrix whose first column is
# y(1:L) and the last row is given by y(L:n)
Y <- trajectory(y, l, k)
# Step 2 (singular value decomposition)
SVD <- svd(Y)
S <- SVD$d
U <- SVD$u
V <- SVD$v
Z <- t(Y)%*%U
# Take II (Reconstruction)
# Step 3 (grouping)
I <- seq(1:m)
Yi <- array(0,c(l, k, m))
YI <- array(0,c(l, k))
for(i in 1:m) {
Yi[, , i] <- U[, i]%*%t(U[, i])%*%Y
YI <- YI + Yi[, , i]
}
# Step 4 (diagonal averaging)
y_ <- Dbar(YI)
answer <- list(fitted.values = y_, PC = Z, trajectory = Y,
series = y, resultant = YI, resultants = Yi,
U = U, S = S, V = V, observations = n, k = k, m = m)
return(answer)
}
```

where Diagonal Averaging (`Dbar`

) is conducted with the auxiliary function

```
Dbar <- function(Y) {
l <- nrow(Y)
k <- ncol(Y)
n <- l + (k - 1)
sum <- rep(0, n)
count <- rep(0, n)
for(L in 2:(k + l)) {
for(i in 1:l) {
for(j in 1:k) {
if(i + j == L) {
count[L - 1] <- count[L - 1] + 1
sum[L - 1] = sum[L - 1] + Y[i, j]
}
}
}
}
ans <- sum / count
}
```

and with the trajectory matrix being constructed using

```
trajectory <- function(x, nrow = length(x)%/%2, ncol = length(x)%/%2) {
Z <- matrix(x[1:ncol], nrow = nrow, ncol = ncol, byrow = T)
for (i in 1:(nrow - 1)) {
Z[i + 1, ] <- x[-(1:i)][1:ncol]}
Z
}
```

The US GDP data used by de Carvalho et al. (2012) can be downloaded from here. We are now ready to conduct the analysis.

```
# Load data and convert it into a time series object
GDP <- read.csv2("us_gdp.csv", header = F)
gdp <- ts(log(GDP), frequency = 4, start = c(1950, 1))
n <- length(gdp)
# Disentangle gdp data into components using ssa
pc <- matrix(0, nrow = n, ncol = 12)
l <- 32
m <- 12
fit <- ssa(gdp, l, m)
for(i in 1:m)
pc[, i] <- Dbar(fit$resultants[, , i])
PC <- list()
for(i in 1:m)
PC[[i]] <- ts(pc[, i], frequency = 4, start = c(1950,1))
```

Our business cycle indicator is given by removing the components which are related with trend (`pc1`

, `pc2`

) and with noise (`pc10`

, `pc11`

, `pc12`

):

`cycle <- PC[[3]] + PC[[4]] + PC[[5]] + PC[[6]] + PC[[7]] + PC[[8]] + PC[[9]]`

The next lines of code yield a comparison of our filter with that of Christiano–Fitzgerald, which can be readily computed using `cffilter`

from the package `mFilter`

:

`cfcycle <- cffilter(gdp, pl = 6, pu = 32, root = T)$cycle`

We now plot the resulting outcomes. We start with the command `rect`

so to draw a sequence of rectangles representing contraction periods as dated by the NBER (National Bureau of Economics Research)

```
# (m, M) - settings which delimit the height of shaded areas
m <- min(c(cycle, cfcycle)) ; M <- max(c(cycle, cfcycle))
plot(cycle, ylim = c(m, M),
main = "Singular Spectrum Filter vs. Christiano-Fitzgerald Filter",
xlab = "Time (in years)", ylab = "")
# calendar definition: S and E respectively denote "start" and "end"
# (certainly not efficient, but helps clarifying what the rectangles below mean)
Sjan = 0/12; Sfev = 1/12; Smar = 2/12; Sapr = 3/12; Smay = 4/12; Sjun = 5/12
Sjul = 6/12; Saug = 7/12; Ssep = 8/12; Soct = 9/12; Snov = 10/12; Sdec = 11/12
Ejan = 1/12; Efev = 2/12; Emar = 3/12; Eapr = 4/12; Emay = 5/12; Ejun = 6/12
Ejul = 7/12; Eaug = 8/12; Esep = 9/12; Eoct = 10/12; Enov = 11/12; Edec = 12/12
rect((Sjul + 1953), m, Emay + 1954, M, col = "gray", border = NA)
rect((Saug + 1957), m, Eapr + 1958, M, col = "gray", border = NA)
rect((Sapr + 1960), m, Efev + 1961, M, col = "gray", border = NA)
rect((Sdec + 1969), m, Enov + 1970, M, col = "gray", border = NA)
rect((Snov + 1973), m, Emar + 1975, M, col = "gray", border = NA)
rect((Sjan + 1980), m, Ejul + 1980, M, col = "gray", border = NA)
rect((Sjul + 1981), m, Enov + 1982, M, col = "gray", border = NA)
rect((Sjul + 1990), m, Emar + 1991, M, col = "gray", border = NA)
rect((Smar + 2001), m, Enov + 2001, M, col = "gray", border = NA)
rect((Sdec + 2007), m, Ejun + 2009, M, col = "gray", border = NA)
lines(cycle)
lines(cfcycle, lty = 2)
```

You may want to add a 'P' and a 'T'—respectively denoting Peak and Trough—at the beginning and at the of the gray rectangles, and this can be done easily using the `R`

command `text`

.

de Carvalho, M., Rodrigues, P. C., and Rua, A. (2012), “Tracking the US Business Cycle with a Singular Spectrum Analysis,” *Economics Letters*, 114, 32–35.