Forecasting GDP with R and dataseries.org


The website dataseries.org aims to be Switzerland’s FRED - a free comprehensive database of Swiss time series. Powered by R and written in Shiny (also using a bit of JavaScript) it allows you to quickly search and explore a large number of data series.

www.dataseries.org

Switzerland’s time series in one place

Similarly to the United States, public data in Switzerland is produced by a large number of different offices, which makes it hard to find any particular series. dataseries.org provides a structured and automatically updated collection of most of these data. We are still working on the data input, but are pretty much complete in the field of Economics.

You can download the data as spreadsheets or graphs, or embed interactive widgets in your website. Alternatively, you can import the data directly into R, using the dataseries package. Install the package from CRAN:

install.packages("dataseries")

and run the ds function with the id argument that you find on the website:

plot(dataseries::ds("GDP.PBRTT.A.R", "ts"), 
     ylab = "mio CHF, at 2010 prices, s. adj.", main = "Gross Domestic Product")

This will give you an R plot of Switzerland’s GDP. (The data is cached, so calling the function again will not re-download until you restart the R session.)


Live Import of Series to R

In the following, I will use data from dataseries.org to produce a live forecast of Switzerland’s GDP. Each day the model is run, it will be ensured that the latest data is used. That way it is possible to produce a transparent and up-to-date forecast. For the following exercise, I will only use tools from R base, but it is of course possible to use the same data in a more advanced modeling framework.

In order to produce a reasonable forecast, we want to track early information on the business cycle, which is mostly survey data. We will use a question on from the SECO Consumer Confidence Survey on current economic performance, the Credit Suisse / Procure Purchasing Managers’ Index and the ETHZ KOF Barometer.

Transforming the data

Getting these indicators from dataseries.org directly into R is easy. Because these data are measured at different frequencies, we need to convert them to the same quarterly frequency as GDP. There are many packages that offer functions for that (e.g., our tempdisagg package has functions to move both to higher or lower frequencies), but I will stick to basic R here:

# Aggregating months to quarters (post updated on May 6, 17)
to_quarterly <- function(x){
  aggregate(x, nfrequency = 4, FUN = mean)
}

pmi <- to_quarterly(dataseries::ds("PMI.SA.PM", "ts"))
kof <- to_quarterly(dataseries::ds("KOF.KFBR", "ts"))
csent <- dataseries::ds("CCI.GEPC", "ts")

A plot of these series shows the common trend in these variables and gives you an indication of the business cycle, which may have turned upward in recent months.

plot(cbind(pmi, kof, csent), main = "Business Cycle Indicators")

Since these series are stationary, our left hand side variable should be stationary as well. This is accomplished by calculating percentage change rates of GDP:

gdp.level <- dataseries::ds("GDP.PBRTT.A.R", "ts")
gdp <- (gdp.level / lag(gdp.level, -1)) - 1   

ARIMA modelling

R’s workhorse for time series modeling is the arima function, which allows you to construct a univariate or multivariate model of GDP growth. Since the data is seasonally adjusted, a simple autoregressive process (AR1) offers a good benchmark:

# AR1
m0 <- arima(gdp, order = c(1, 0, 0))
fct0 <- predict(m0, n.ahead = 1)$pred
# GDP Growth Q1: +0.3 

If you need advice on which ARIMA model to choose, the information criterions, accessed by the R functions AIC or BIC, can help you to choose a model. Simply take the model with lowest information criterion. The auto.arima function from the forecast package also allows you to do the selection automatically.

We can include our series individually or jointly and estimate a range of different models. A good model (in terms of the AIC information criterion) is the following, which uses PMI and KOF data (but not consumer sentiment data):

# PMI, KOF
dta <- window(cbind(pmi, kof), start = start(pmi), end = end(pmi))
m1 <- arima(window(gdp, start = start(dta)), 
            xreg = window(dta, end = end(gdp)))
fct1 <- predict(m1, n.ahead = 1, 
                newxreg = window(dta, start = tsp(gdp)[2] + 0.25,
                                 end = tsp(gdp)[2] + 0.25)
                )$pred
# GDP Growth Q1: +0.7 

The model’s forecast for the first quarter of 2017 is 0.7 - a value that hasn’t been reached for more than two years.

A factor model

If you have multiple indicators at hand, a common problem is multicollinearity, the fact that indicators are correlated, and therefore too many indicators deteriorate the quality of the model estimation.

An easy fix is to use a factor model, where the indicators are summarized in a few factors, which can be calculated by principal components (see Stock and Watson 2002):

# PMI, KOF, Consumer Sentiment, first Principal Component
pca <- prcomp(window(cbind(pmi, kof, csent), start = start(pmi), 
                     end = tsp(gdp)[2] + 0.25),
              scale. = TRUE)
dta.pca <- ts(pca$x[, 'PC1'], start = start(pmi), frequency = 4)

m2 <- arima(window(gdp, start = start(dta)), 
            xreg = window(dta.pca, end = end(gdp)))
fct2 <- predict(m2, n.ahead = 1, 
                newxreg = window(dta.pca, start = tsp(gdp)[2] + 0.25)
                )$pred
# GDP Growth Q1: +0.7 

Again, we get a forecast value of 0.7. Overall, survey data indicates that the economy is well on track. Let’s do a graphical comparison of our forecasts:

# skeletons to include forecasts 
gdp.fct0  <- window(gdp, extend = TRUE, end = tsp(gdp)[2] + 0.25)
gdp.fct1 <- gdp.fct2 <- gdp.fct0

# plug forecasts into skeletons
window(gdp.fct0, start = end(gdp.fct0)) <- fct0
window(gdp.fct1, start = end(gdp.fct1)) <- fct1
window(gdp.fct2, start = end(gdp.fct2)) <- fct2

ts.plot(window(cbind(gdp, gdp.fct0, gdp.fct1, gdp.fct2), start = 2010), 
        col = 1:4, ylab = "quarterly growth rates, s. adj.", 
        main = "GDP Forecasts")
legend("topright", 
       legend = c("GDP Growth Rate", "AR 1 Forecast", 
                  "PMI, KOF", "Principal Component"), 
       lty = 1, col = 1:4, bty = "n")

Publication of first quarter GDP is on June 1, 2017. See you in a month!

seasonal 1.3: A Better Way to Seasonal Adjustment Diagnostics


The R package seasonal makes it easy to use X-13ARIMA-SEATS, the seasonal adjustment software by the United States Census Bureau. Thanks to the x13binary package, installing it from CRAN is now as easy as installing any other R package:

install.packages("seasonal")

The latest version 1.3 comes with a new udg function and a customizable summary method, which give power users of X-13 a convenient way to check the statistics that are of their interest. For a full list of changes, see the package NEWS.

A generalized way to access diagnostics

Version 1.3 offers a generalized way to access diagnostic statistics. In seasonal, it was always possible to use all options of X-13 and access all output series. Now it is easy to access all diagnostics as well.

The main new function is udg, named after the X-13 output file which it is reading. Let us start with a simple call to seas (the main function of the seasonal package) that uses the X-11 seasonal adjustment method:

m <- seas(AirPassengers, x11 = "")
udg(m)

The udg function returns a list containing 357 named diagnostics. They are properly type-converted, so they can be directly used for further analysis within R.

If we ask for a specific statistic, such as the popular X-11 M statistics, the result will be simplified to a numeric vector (see ?udg for additional options):

udg(m, c("f3.m01", "f3.m02", "f3.m03", "f3.m04"))

## f3.m01 f3.m02 f3.m03 f3.m04 
##  0.041  0.042  0.000  0.283

There are also some new wrappers for commonly used statistics, such as AIC, BIC, logLik or qs, which use the udg function.

A customizable summary

The new functionality paves the way for a customizable summary method for seas objects. For example, if we want to add the M statistics for X-11 adjustments to the summary, we can write:

summary(m, stats = c("f3.m01", "f3.m02", "f3.m03", "f3.m04"))

## Call:
## seas(x = AirPassengers, x11 = "")
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## Weekday           -0.0029497  0.0005232  -5.638 1.72e-08 ***
## Easter[1]          0.0177674  0.0071580   2.482   0.0131 *  
## AO1951.May         0.1001558  0.0204387   4.900 9.57e-07 ***
## MA-Nonseasonal-01  0.1156204  0.0858588   1.347   0.1781    
## MA-Seasonal-12     0.4973600  0.0774677   6.420 1.36e-10 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## X11 adj.  ARIMA: (0 1 1)(0 1 1)  Obs.: 144  Transform: log
## AICc: 947.3, BIC: 963.9  QS (no seasonality in final):    0  
## Box-Ljung (no autocorr.): 26.65   Shapiro (normality): 0.9908  
## f3.m01: 0.041  f3.m02: 0.042  f3.m03: 0  f3.m04: 0.283 

Note the new line at the end, which shows the M statistics.

If we want to routinely consider these statistics in our summary, we can set the seas.stats option accordingly:

options(seas.stats = c("f3.m01", "f3.m02", "f3.m03", "f3.m04"))

This will change the default behavior, and

summary(m)

will return the same output as above. To restore the default behavior, set the option back to NULL.

options(seas.stats = NULL)

Like Peanut Butter and Jelly: x13binary and seasonal


This post was written by Dirk Eddelbuettel and Christoph Sax and posted by both author’s respective blogs.

The seasonal package by Christoph Sax brings a very featureful and expressive interface for working with seasonal data to the R environment. It uses the standard tool of the trade: X-13ARIMA-SEATS. This powerful program is provided by the statisticians of the US Census Bureau based on their earlier work (named X-11 and X-12-ARIMA) as well as the TRAMO/SEATS program by the Bank of Spain. X-13ARIMA-SEATS is probably the best known tool for de-seasonalization of timeseries, and used by statistical offices around the world.

Flattening the Learning Curve

Sadly, it also has a steep learning curve. One interacts with a basic command-line tool which users have to download, install and properly reference (by environment variables or related means). Each model specification has to be prepared in a special ‘spec’ file that uses its own, cumbersome syntax.

As seasonal provides all the required functionality to use X-13ARIMA-SEATS from R — see the very nice seasonal demo site — it still required the user to manually deal with the X-13ARIMA-SEATS installation.

So we decided to do something about this. A pair of GitHub repositories provide both the underlying binary in a per-operating system form (see x13prebuilt) as well as a ready-to- use R package (see x13binary) which uses the former to provide binaries for R. And the latter is now on CRAN as package x13binary ready to be used on Windows, OS-X or Linux.

Automatic Deployment of X-13

And the seasonal package (in version 1.2.0 – now on CRAN – or later) automatically makes use of it. Installing seasaonal and x13binary in R is now as easy as:

install.packages("seasonal")

which opens the door for effortless deployment of powerful deasonalization. By default, the principal function of the package employs a number of automated techniques that work well in most circumstances. For example, the following code produces a seasonal adjustment of the latest data of US retail sales (by the Census Bureau) downloaded from Quandl:

library(seasonal) 

url <- "https://www.quandl.com/api/v3/datasets/USCENSUS/BI_MARTS_44000_SM.csv?order=asc"
rs <- ts(read.csv(url)$Value/1e3, start = c(1992, 1), frequency = 12)

m1 <- seas(rs)

plot(m1, main = "Retail Trade: U.S. Total Sales", ylab = "USD (in Billions)")

This tests for log-transformation, performs an automated ARIMA model search, applies outlier detection, tests and adjusts for trading day and easter effects, and invokes the SEATS method to perform seasonal adjustment. And this is how the adjusted series looks like:

US Retail Sales, seasonally adjusted

Fig 1: US Retail Sales, seasonally adjusted

Of course, you can access all available options of X-13ARIMA-SEATS as well. Here is an example where we adjust the latest data for Chinese exports (as tallied by the US FED), taking into account the different effects of Chinese New Year before, during and after the holiday:

url <- "https://www.quandl.com/api/v3/datasets/FRED/VALEXPCNM052N.csv?order=asc"
xp <- ts(read.csv(url)$VALUE/1e9, start = c(1981, 1), frequency = 12)

m2 <- seas(window(xp, start = 2000),
          xreg = cbind(genhol(cny, start = -7, end = -1, center = "calendar"), 
                       genhol(cny, start = 0, end = 7, center = "calendar"), 
                       genhol(cny, start = 8, end = 21, center = "calendar")
          ),
          regression.aictest = c("td", "user"),
          regression.usertype = "holiday")

plot(m2, main = "Goods, Value of Exports for China", ylab = "USD (in Billions)")

which generates the following chart demonstrating a recent flattening in export activity measured in USD.

Chinese Exports, seasonally adjusted

Fig 2: Chinese Exports, seasonally adjusted

We hope this simple examples illustrates both how powerful a tool X-13ARIMA-SEATS is, but also just how easy it is to use X-13ARIMA-SEATS from R now that we provide the x13binary package automating its installation.

Christoph Sax Data Analytics LLC
Langstrasse 191, 8005 Zurich, Switzerland
Phone: +41 (0)43 540 26 91
info@christophsax.com

© 2016 Christoph Sax Data Analytics