Bank of England Fan Charts in R

[This article was first published on Guy Abel » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

I managed to catch David Spiegelhalter’s Tails You Win on BBC iplayer last week. I missed it the first time round, only for my parents on my last visit home to tell me about a Statistician jumping out of a plane on TV. It was a great watch. Towards the end I spotted some fan charts used by the Bank of England to illustrate uncertainty in their forecasts, similar to this one:

Bank of England February 2013 CPI Fan Chart.
Obtained from http://www.bankofengland.co.uk/publications/Pages/inflationreport/irfanch.aspx

They discussed how even in the tails of their GDP predictive distribution they missed the financial crisis by a long shot. This got me googling, trying and find you how they made the plots, something that (also) completely passed me by when I put together my fanplot package for R. As far as I could tell they did them in Excel, although (appropriately) I am not completely certain. There are also MATLAB files that can create fan charts. Anyhow, I thought I would have a go at replicating a Bank of England fan chart in R….

Split Normal (Two-Piece) Normal Distribution.

The Bank of England produce fan charts of forecasts for CPI and GDP in their quarterly Inflation Reports. They also provide data for each report, in the form of mean, uncertainty and a skewness indicator. These three parameters are from a split-normal distribution (the BOE call it a Two-Piece Normal Distribution):
spn
where
spnsig
so if you know \sigma and \gamma you can derive \sigma_1 and \sigma_2. As no split normal distribution existed in R, I added routines for a density, distribution and quantile function, plus a random generator, to a new version (2.1) of the fanplot package. I used the formula in Julio (2007) to code each of the three functions, and checked the results against those from the fan chart MATLAB code.

Fan Chart Plots for CPI.

Once I had the qsplitnorm function working properly, producing the fan chart plot in R was pretty straight-forward. First I downloaded the past CPI data (percentage change over 12 months) from ONS:

> cpi <- read.csv("C:/Users/…/cpi.csv")
> cpi
   Year Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1  1997 2.1 1.9 1.7 1.6 1.6 1.7 2.0 2.0 1.8 1.9 1.9 1.7
2  1998 1.5 1.6 1.7 1.8 2.0 1.7 1.4 1.3 1.4 1.4 1.4 1.6
3  1999 1.6 1.4 1.7 1.5 1.3 1.3 1.3 1.2 1.2 1.1 1.2 1.1
4  2000 0.8 0.9 0.6 0.6 0.5 0.8 0.9 0.6 1.0 1.0 1.1 0.8
5  2001 0.9 0.8 0.9 1.2 1.7 1.7 1.4 1.8 1.3 1.2 0.8 1.1
6  2002 1.6 1.5 1.5 1.4 0.8 0.6 1.1 1.0 1.0 1.4 1.5 1.7
7  2003 1.3 1.6 1.5 1.4 1.3 1.1 1.3 1.4 1.4 1.4 1.3 1.3
8  2004 1.4 1.3 1.1 1.1 1.5 1.6 1.4 1.3 1.1 1.2 1.5 1.7
9  2005 1.6 1.7 1.9 1.9 1.9 2.0 2.3 2.4 2.5 2.3 2.1 1.9
10 2006 1.9 2.0 1.8 2.0 2.2 2.5 2.4 2.5 2.4 2.4 2.7 3.0
11 2007 2.7 2.8 3.1 2.8 2.5 2.4 1.9 1.8 1.8 2.1 2.1 2.1
12 2008 2.2 2.5 2.5 3.0 3.3 3.8 4.4 4.7 5.2 4.5 4.1 3.1
13 2009 3.0 3.2 2.9 2.3 2.2 1.8 1.8 1.6 1.1 1.5 1.9 2.9
14 2010 3.5 3.0 3.4 3.7 3.4 3.2 3.1 3.1 3.1 3.2 3.3 3.7
15 2011 4.0 4.4 4.0 4.5 4.5 4.2 4.4 4.5 5.2 5.0 4.8 4.2
16 2012 3.6 3.4 3.5 3.0 2.8 2.4 2.6 2.5 2.2 2.7 2.7 2.7
> 
> # extract quarterly data in Feb, Apr, Aug and Nov
> mn <- month.abb[seq(2, 11, 3)]
> cpi <- cpi[, mn]
> # covert to ts
> cpi <- ts(c(t(cpi)), start = 1997, frequency = 4)
> cpi
     Qtr1 Qtr2 Qtr3 Qtr4
1997  1.9  1.6  2.0  1.9
1998  1.6  2.0  1.3  1.4
1999  1.4  1.3  1.2  1.2
2000  0.9  0.5  0.6  1.1
2001  0.8  1.7  1.8  0.8
2002  1.5  0.8  1.0  1.5
2003  1.6  1.3  1.4  1.3
2004  1.3  1.5  1.3  1.5
2005  1.7  1.9  2.4  2.1
2006  2.0  2.2  2.5  2.7
2007  2.8  2.5  1.8  2.1
2008  2.5  3.3  4.7  4.1
2009  3.2  2.2  1.6  1.9
2010  3.0  3.4  3.1  3.3
2011  4.4  4.5  4.5  4.8
2012  3.4  2.8  2.5  2.7

Then I read in the parameters for the forecast distribution from the Bank of England:

> cpiboe <- read.csv("C:/Users/…/cpiboe.csv")
> cpiboe
   year quarter mode median mean uncertainty skewness
1  2013      Q1 2.73   2.73 2.73        0.61        0
2  2013      Q2 2.92   2.92 2.92        0.88        0
3  2013      Q3 3.22   3.22 3.22        1.11        0
4  2013      Q4 3.13   3.13 3.13        1.27        0
5  2014      Q1 2.95   2.95 2.95        1.34        0
6  2014      Q2 2.82   2.82 2.82        1.37        0
7  2014      Q3 2.53   2.53 2.53        1.39        0
8  2014      Q4 2.41   2.41 2.41        1.42        0
9  2015      Q1 2.32   2.32 2.32        1.48        0
10 2015      Q2 2.23   2.23 2.23        1.49        0
11 2015      Q3 2.13   2.13 2.13        1.50        0
12 2015      Q4 2.01   2.01 2.01        1.52        0
13 2016      Q1 1.96   1.96 1.96        1.52        0

Using the qsplitnorm function I then derived the CPI values for a sequence of probabilities between 0 and 1:

> library("fanplot")
> x <- 1:99/100
> k <- nrow(cpiboe)
> xx <- matrix(NA, nrow = 99, ncol = k)
> for (i in 1:k) {
    xx[, i] <- qsplitnorm(x, mean = cpiboe$mean[i],
                             sd = (cpiboe$uncertainty[i])^2,
                             skew = cpiboe$skewness[i])
  }
> head(xx)
         [,1]     [,2]      [,3]        [,4]         [,5]        [,6]       [,7]       [,8]      [,9]     [,10]     [,11]     [,12]     [,13]
[1,] 1.864366 1.118476 0.3537068 -0.62216649 -1.227190243 -1.54632232 -1.9647367 -2.2808479 -2.775632 -2.934725 -3.104283 -3.364794 -3.414794
[2,] 1.965800 1.329577 0.6895760 -0.18249162 -0.737711544 -1.03468133 -1.4380483 -1.7311793 -2.178532 -2.329528 -2.490935 -2.734981 -2.784981
[3,] 2.030157 1.463513 0.9026742  0.09646799 -0.427153003 -0.71006152 -1.1038813 -1.3824322 -1.799690 -1.945550 -2.101786 -2.335386 -2.385386
[4,] 2.078570 1.564269 1.0629797  0.30631844 -0.193531910 -0.46586269 -0.8525006 -1.1200834 -1.514703 -1.656698 -1.809044 -2.034785 -2.084785
[5,] 2.117950 1.646225 1.1933758  0.47701559 -0.003499173 -0.26722577 -0.6480217 -0.9066829 -1.282887 -1.421740 -1.570921 -1.790270 -1.840270
[6,] 2.151469 1.715983 1.3043635  0.62230567  0.158248534 -0.09815456 -0.4739781 -0.7250455 -1.085576 -1.221753 -1.368241 -1.582149 -1.632149

These values are ready to be passed to the pn function to calculate the percentiles (in a pn.ts type object) for plotting.

> # save time of last boe data point
> boe0 <- tsp(cpi)[2]
> # create pn.ts type object ready for fan function
> xx.pn <- pn(xx, start = boe0, frequency = 4, anchor = tail(cpi, 1))
> # view head of object for 5th, 10th, 50th, 90th and 95th percentiles
> head(xx.pn, p = c(5, 10, 50))
            5%       10%  50%      90%      95%
[1,] 2.7000000 2.7000000 2.70 2.700000 2.700000
[2,] 2.1481169 2.2695140 2.73 3.190486 3.311883
[3,] 1.7090075 1.9616546 2.92 3.878345 4.130992
[4,] 1.2932647 1.6952358 3.22 4.744764 5.146735
[5,] 0.6077767 1.1339833 3.13 5.126017 5.652223
[6,] 0.1420738 0.7278861 2.95 5.172114 5.757926
> # plot the data
> plot(cpi, type = "l", xlim = c(1997, 2016), ylim = c(-2, 7), las = 1)
> # add fan
> fan(xx.pn)

Note, setting the anchor to the last data point joins the fan to the past data line.

This plot is based on the default arguments (colours, percentiles, lines and labels) in the pn and fan functions. The Bank of England fan charts are somewhat different. They use a courser set of colours and do not show the tails of the distribution. These can be derived by specifying which percentiles to estimate from the pn function (by default, as above, percentiles 1 to 100 are calculated for each time period):

> # percentiles in steps of 5, without tails.
> xx.pn <- pn(xx, p = seq(15, 85, 5), start = boe0, frequency = 4, anchor = tail(cpi, 1))
> head(xx.pn)
          15%      20%      25%      30%      35%      40%      45%  50%      55%      60%      65%      70%      75%      80%      85%
[1,] 2.700000 2.700000 2.700000 2.700000 2.700000 2.700000 2.700000 2.70 2.700000 2.700000 2.700000 2.700000 2.700000 2.700000 2.700000
[2,] 2.355276 2.424691 2.484817 2.539120 2.589621 2.637650 2.684180 2.73 2.775820 2.822350 2.870379 2.920880 2.975183 3.035309 3.104724
[3,] 2.140140 2.284604 2.409734 2.522748 2.627848 2.727804 2.824641 2.92 3.015359 3.112196 3.212152 3.317252 3.430266 3.555396 3.699860
[4,] 1.979213 2.209060 2.408148 2.587957 2.755176 2.914209 3.068281 3.22 3.371719 3.525791 3.684824 3.852043 4.031852 4.230940 4.460787
[5,] 1.505728 1.806614 2.067232 2.302614 2.521514 2.729700 2.931390 3.13 3.328610 3.530300 3.738486 3.957386 4.192768 4.453386 4.754272
[6,] 1.141740 1.476708 1.766848 2.028892 2.272588 2.504356 2.728892 2.95 3.171108 3.395644 3.627412 3.871108 4.133152 4.423292 4.758260

The xx.pn object can then be plotted, having set up a matching colour scheme. I also add shading for the forecast area, alternative axis and lines for the target rate and two year ahead period to match the Bank of England plot:

> # create colour palette
> pal <- colorRampPalette(c("tomato", "white"))
> fancol <- pal(ncol(xx.pn)/2)
> # plot data
> plot(cpi, type = "l", lwd = 3, col = "tomato",
            xlim = c(2008, 2016), ylim = c(-2, 7),
            las = 1, xaxt = "n", yaxt = "n")
> # add shading for forecast area
> lim <- par("usr")
> rect(boe0, lim[3] - 1, 2017, lim[4] + 1, border = "gray90", col = "gray90")
> # add fan
> fan(xx.pn, fan.col = fancol, txt = NA, ln.col = NA)
> # add axis
> axis(2, at = -2:7, las = 2, tcl = 0.5, labels = FALSE)
> axis(4, at = -2:7, las = 2, tcl = 0.5)
> axis(1, at = 2008:2016, tcl = 0.5)
> axis(1, at = seq(2008, 2016, 0.25), labels = FALSE, tcl = 0.2)
> # add target line
> abline(h = 2)
> # add two year line
> abline(v = boe0 + 2, lty = 2)


To leave a comment for the author, please follow the link and comment on their blog: Guy Abel » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)