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.

Singular Spectrum Business Cycle Analysis

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.


References

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.