Bank of England Fan Charts in R
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
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):
where
so if you know 
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)
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.