Meeting the Sustainable Development Goals leads to lower world population growth

A couple of months ago, a paper I worked on with co-authors from the Wittgenstein Centre was published in the Proceedings for the National Academy of Science of the United States of America. The significance blurb and abstract are below.

Abel, G.J., Bakarat, B., K.C., S. and Lutz, W. (2016). Meeting the Sustainable Development Goals leads to lower world population growth. Proceedings of the National Academy of Sciences. 113 (50)  14294–14299

Continue reading “Meeting the Sustainable Development Goals leads to lower world population growth”

Advertisements

Updated Circular Plots for Directional Bilateral Migration Data

I have had a few emails recently regarding plots from my new working paper on global migration flows, which has received some media coverage here, here and here. The plots were created using Zuguang Gu’s excellent circlize package and are a modified version of those discussed in an earlier blog post. In particular, I have made four changes:

  1. I have added arrow heads to better indicate the direction of flows, following the example in Scientific American.
  2. I have reorganized the sectors on the outside of the circle so that in each the outflows are plotted first (largest to smallest) followed by the inflows (again, in size order). I prefer this new layout (previously the inflows were plotted first) as it allows the time sequencing of migration events (a migrant has to leave before they can arrive) to match up with the natural tendency for most to read from left to right.
  3. I have cut out the white spaces that detached the chords from the outer sector. To my eye, this alteration helps indicate the direction of the flow and gives a cleaner look.
  4. I have kept the smallest flows in the plot, but plotted their chords last, so that the focus is maintained on the largest flows. Previously smaller flows were dropped according to an arbitrary cut off, which meant that the sector pieces on the outside of the circle no longer represented the total of the inflows and outflows.

Combined, these four modifications have helped me when presenting the results at recent conferences, reducing the time I need to spend explaining the plots and avoiding some of the confusion that occasionally occurred with the direction of the migration flows.

If you would like to replicate one of these plot, you can do so using estimates of the minimum migrant transition flows for the 2010-15 period and the demo R script in my migest package;

# install.packages("migest")
# install.packages("circlize")
library("migest")
demo(cfplot_reg2, package = "migest", ask = FALSE)

which will give the following output:

Estimated Global Migration Flows 2010-15

The code in the demo script uses the chordDiagram function, based on a recent update to the circlize package (0.3.7). Most likely you will need to either update or install the package (uncomment the install.packages lines in the code above).

If you want to view the R script in detail to see which arguments I used, then take a look at the demo file on GitHub here. I provide some comments (in the script, below the function) to explain each of the argument values.

Save and view a PDF version of the plot (which looks much better than what comes up in my non-square RStudio plot pane) using:

dev.copy2pdf(file ="cfplot_reg2.pdf", height=10, width=10)
file.show("cfplot_reg2.pdf")

Shiny App for the Wittgenstein Centre Population Projections

A few weeks ago a new version of the the Wittgenstein Centre Data Explorer was launched. The data explorer is intended to disseminate the results of a recent global population projection exercise which uniquely incorporates level of education (as well as age and sex) and the scientific input of more than 500 population experts around the world. Included are the projected populations used in the 5th assessment report of the Intergovernmental Panel on Climate Change (IPCC).

Untitled

Over the past year or so I have been working (on and off) with the data lab team to create a shiny app, on which the data explorer is based. All the code and data is available on my github page. Below are notes to summarise some of the lessons I learnt:

1. Large data

We had a pretty large amount of data to display (31 indicators based on up to 7 scenarios x 26 time periods x 223 geographical areas x 21 age groups x 2 genders x 7 education categories)… so somewhere over 8 million rows for some indicators. Further complexity was added by the fact that some indicators were by definition not available for some dimensions of the data, for example, population median age is not available by age group. The size and complexity meant that data manipulations were a big issue. Using read.csv to load the data didn’t really cut the mustard, taking over 2 minutes when running on the server. The fantastic saves package and ultra.fast=TRUE argument in the loads function came to the rescue, alongside some pre-formatting to avoid as much joining and reshaping of the data on the server as possible. This cut load times to a couple of seconds at most, and allowed the app to work with the indicator variables on the fly as demanded by the user selections. Once the data was in, the more than awesome dplyr functions finished the data manipulations jobs in style. I am sure there is some smarter way to get everything running a little bit quicker than it does now, but I am pretty happy with the present speed, given the initial waiting times.

2. googleVis and gvisMerge

It’s a demographic data explorer, which means population pyramids have to pop-up somewhere. We needed pyramids that illustrate population sizes by education level, on top of the standard age and sex breakdown. Static versions of the education pyramids in the explorer have previously been used by my colleagues to illustrate past and future populations. For the graphic explorer I created some interactive versions, for side-by-side comparisons over time and between countries, and which also have some tool tip features. These took a little while to develop. I played with ggvis but couldn’t get my bar charts to go horizontal. I also took a look at some other functions for interactive pyramids but I couldn’t figure out a way to overlay the educational dimension. I found a solution by creating gender specific stacked bar charts from gvisBarChart in the googleVis package and then gvisMerge to bring them together in one plot. As with the data tables, they take a second or so render, so I added a withProgress bar to try and keep the user entertained.

I could not figure out a way in R to convert the HTML code outputted by the gvisMerge function to a familiar file format for users to download. Instead I used a system call to the wkhtmltopdf program to return PDF or PNG files. By default, wkhtmltopdf was a bit hit and miss, especially with converting the more complex plots in the maps to PNG files. I found setting --enable-javascript --javascript-delay 2000 helped in many cases.

3. The shiny user community

I asked questions using the shiny tag on stackoverflow and the shiny google group a number of times. A big thank you to everyone who helped me out. Browsing through other questions and answers was also super helpful. I found this question on organising large shiny code particularly useful. Making small changes during the reviewing process became a lot easier once I broke the code up across multiple .R files with sensible names.

4. Navbar Pages

When I started building the shiny app I was using a single layout with a sidebar and tabbed pages to display data and graphics (using tabsetPanel()), adding extra tabs as we developed new features (data selection, an assumption data base, population pyramids, plots of population size, maps, FAQ’s, etc, etc.). As these grew, the switch to the new Navbar layout helped clean up the appearance and provide a better user experience, where you can move between data, graphics and background information using the bar at the top of page.

5. Shading and link buttons

I added some shading and buttons to help navigate through the data selection and between different tabs. For the shading I used cssmatic.com to generate the colour of a fluidRow background. The code generated there was copy and pasted into a tags$style element for my defined row myRow1, as such;

library(shiny)
runApp(list(
  ui = shinyUI(fluidPage(
    br(),
    fluidRow(
      class = "myRow1", 
      br(),
      selectInput('variable', 'Variable', names(iris))
    ),
    tags$style(".myRow1{background: rgba(212,228,239,1); 
                background: -moz-linear-gradient(left, rgba(212,228,239,1) 0%, rgba(44,146,208,1) 100%);
                background: -webkit-gradient(left top, right top, color-stop(0%, rgba(212,228,239,1)), color-stop(100%, rgba(44,146,208,1)));
                background: -webkit-linear-gradient(left, rgba(212,228,239,1) 0%, rgba(44,146,208,1) 100%);
                background: -o-linear-gradient(left, rgba(212,228,239,1) 0%, rgba(44,146,208,1) 100%);
                background: -ms-linear-gradient(left, rgba(212,228,239,1) 0%, rgba(44,146,208,1) 100%);
                background: linear-gradient(to right, rgba(212,228,239,1) 0%, rgba(44,146,208,1) 100%);
                filter: progid:DXImageTransform.Microsoft.gradient( startColorstr='#d4e4ef', endColorstr='#2c92d0', GradientType=1 );
                border-radius: 10px 10px 10px 10px;
                -moz-border-radius: 10px 10px 10px 10px;
                -webkit-border-radius: 10px 10px 10px 10px;}")
  )),
  server = function(input, output) {
  }
))

I added some buttons to help novice users switch between tabs once they had selected or viewed their data. It was a little tougher to implement than the shading, and in the end I need a little help. I used bootsnipp.com to add some icons and define the style of the navigation buttons (using the tags$style element again).

That is about it for the moment. I might add a few more notes to this post as they occur to me… I would encourage anyone who is tempted to learn shiny to take the plunge. I did not know JavaScript or any other web languages before I started, and I still don’t… which is the great beauty of the shiny package. I started with the RStudio tutorials, which are fantastic. The R code did not get a whole lot more complex than what I learnt there, even though the shiny app is quite large in comparisons to most others I have seen.

Any comments or suggestions for improving website are welcome.

Circular Migration Flow Plots in R

Please see this blog post on updated version of circular plots for migration flows, based on global estimates for 2010-15.

A article of mine was published in Science today. It introduces estimates for bilateral global migration flows between all countries. The underlying methodology is based on the conditional maximisation routine in my Demographic Research paper. However, I tweaked the demographic accounting which ensures the net migration in the estimated migration flow tables matches very closely to the net migration figures from the United Nations.

My co-author, Nikola Sander, developed some circular plots for the paper based on circos in perl. A couple of months back, after the paper was already in the submission process, I figured out how to replicate these plots in R using the circlize package. Zuguang Gu, the circlize package developer was very helpful, responding quickly (and with examples) to my emails.

To demonstrate, I have put two demo files in my migest R package. For the estimates of flows by regions, users can hopefully replicate the plots (so long as the circlize and plyr packages are installed) using:

library("migest")
demo(cfplot_reg, package = "migest", ask = FALSE)

It should result in the following plot:
cfplot_reg
The basic idea of the plot is to show simultaneously the relative size of estimated flows between regions. The origins and destinations of migrants are represented by the circle’s segments, where nearby regions are positioned close to each other. The size of the estimated flow is indicated by the width of the link at its bases and can be read using the tick marks (in millions) on the outside of the circle’s segments. The direction of the flow is encoded both by the origin colour and by the gap between link and circle segment at the destination.

You can save the PDF version of the plot (which looks much better than what comes up in my R graphics device) using:

dev.copy2pdf(file = "cfplot_reg.pdf", height=10, width=10)

If you want to view the R script:

file.show(system.file("demo/cfplot_reg.R", package = "migest"))

In Section 5 of our Vienna Institute of Demography Working Paper I provide a more detailed breakdown for the R code in the demo files.

A similar demo with slight alterations to the labelling is also available for a plot of the largest country to country flows:

demo(cfplot_nat, package = "migest", ask = FALSE)

cfplot_nat

If you are interested in the estimates, you can fully explore in the interactive website (made using d3.js) at http://global-migration.info/. There is also a link on the website to download all the data. Ramon Bauer has a nice blog post explaining the d3 version.

Publication Details:

Abel, G.J. and Sander, N. (2014). Quantifying Global International Migration Flows. Science. 343 (6178), 1520–1522.

Widely available data on the number of people living outside of their country of birth do not adequately capture contemporary intensities and patterns of global migration flows. We present data on bilateral flows between 196 countries from 1990 through 2010 that provide a comprehensive view of international migration flows. Our data suggest a stable intensity of global 5-year migration flows at ~0.6% of world population since 1995. In addition, the results aid the interpretation of trends and patterns of migration flows to and from individual countries by placing them in a regional or global context. We estimate the largest movements to occur between South and West Asia, from Latin to North America, and within Africa.

Forecasting Environmental Immigration to the UK

A couple of months ago, a paper I worked on with co-authors from the Centre of Population Change was published in Population and Environment. It summarised work we did as part of the UK Government Office for Science Foresight project on Migration and Global Environmental Change. Our aim was to build expert based forecasts of environmental immigrants to the UK. We conducted a Delphi survey of nearly 30 migration experts from academia, the civil service and non-governmental organisations to obtain estimates on the future levels of immigration to the UK in 2030 and 2060 with uncertainty. We also asked them what proportion of current and future immigration are/will be environmental migrants. The results were incorporated into a set of model averaged Bayesian time series models through prior distributions on the mean and variance terms.

The plots in the journal article got somewhat butchered during the publication process. Below is the non-butchered version for the future immigration to the UK alongside the past immigration data from the Office of National Statistics.
imm2
At first, I was a bit taken aback with this plot. A few experts thought there were going to be some very high levels of future immigration which cause the rather striking large upper tail. However, at a second glance, the central percentiles show a gentle decrease where these is only (approximately) a 30% chance of an increase in future migration from the 2010 level throughout the forecast period.

The expert based forecast for total immigration was combined with the responses to questions on the proportion of environmental migrants, to obtain an estimate on both the current level of environmental migration (which is not currently measured) and future levels:
env4

As is the way with these things, we came across some problems in our project. The first, was with the definition of an environmental migrant, which is not completely nailed on in the migration literature. As a result the part of the uncertainty in the expert based forecasts are reflective of not only the future level but also of the measure itself. The second was with the elicitation of uncertainty. We used a Likert type scale, which caused some difficulties even during the later round of the Delphi survey. If I was to do over, then this I reckon problem could be much better addressed by getting experts to visualise their forecast fans in an interactive website, perhaps creating a shiny app with the fanplot package. Such an approach would result in smoother fans than those in the plots above, which were based on interpolations from expert answers at only two points of time in the future (2030 and 2060).

Publication Details:

Abel, G.J., Bijak, J., Findlay, A.M., McCollum, D. and Wiśniowski, A. (2013). Forecasting environmental migration to the United Kingdom: An exploration using Bayesian models. Population and Environment. 35 (2), 183–203

Over the next 50 years, the potential impact of environmental change on human livelihoods could be considerable, with one possible consequence being increased levels of human mobility. This paper explores how uncertainty about the level of immigration to the United Kingdom as a consequence of environmental factors elsewhere may be forecast using a methodology involving Bayesian models. The conceptual understanding of forecasting is advanced in three ways. First, the analysis is believed to be the first time that the Bayesian modelling approach has been attempted in relation to environmental mobility. Second, the paper considers the expediency of this approach by comparing the responses to a Delphi survey with conventional expectations about environmental mobility in the research literature. Finally, the values and assumptions of the expert evidence provided in the Delphi survey are interrogated to illustrate the limited set of conditions under which forecasts of environmental mobility, as set out in this paper, are likely to hold.

Global Bilateral International Migration Flows

A few months ago, Demographic Research published my paper on estimating global migration flow tables. In the paper I developed a method to estimate international migrant flows, for which there is limited comparable data, to matches changes in migrant stock data, which are more widely available. The result was bilateral tables of estimated international migrant transitions between 191 countries for four decades, which I believe are a first of kind. The estimates in an excel spreadsheet are available as a additional file on the journal website. The abstract and citation details are at the bottom of this post.

My migest R package contains the ffs function for the flows-from-stock method used in the paper. To demonstrate, consider two hypothetical migrant stock tables I use in the paper, where rows represent place of birth and columns represent place of residence. The first stock table represents the distributions of migrant stocks at the start of the period. The second represents the distributions at the end of the period.

> # create P1 and P2 stock tables
> dn <- LETTERS[1:4]
> P1 <- matrix(c(1000, 100, 10, 0,
+                55, 555, 50, 5, 
+                80, 40, 800, 40, 
+                20, 25, 20, 200), 
+              nrow=4, ncol=4, byrow = TRUE,
+              dimnames = list(pob = dn, por = dn))
> P2 <- matrix(c(950, 100, 60, 0, 
+                80, 505, 75, 5, 
+                90, 30, 800, 40,
+                40, 45, 0, 180),
+              nrow=4, ncol=4, byrow = TRUE,
+              dimnames = list(pob = dn, por = dn))
> # display with row and col totals
> addmargins(P1)
     por
pob      A   B   C   D  Sum
  A   1000 100  10   0 1110
  B     55 555  50   5  665
  C     80  40 800  40  960
  D     20  25  20 200  265
  Sum 1155 720 880 245 3000
> addmargins(P2)
     por
pob      A   B   C   D  Sum
  A    950 100  60   0 1110
  B     80 505  75   5  665
  C     90  30 800  40  960
  D     40  45   0 180  265
  Sum 1160 680 935 225 3000

When estimating flows from stock data, a good demographer should worry about births and deaths over the period as these can have substantial impacts on changes in populations over time. In the simplest example using the above hypothetical example above, I set births and deaths to zero (implied by the equal row totals, the sum of populations by their place of birth) in each stock table. In any case I need to create some vectors to pass this information to the ffs function.

> # no births and deaths
> b <- rep(0, 4)
> d <- rep(0, 4)

We can then pass the stock tables, births and deaths to the ffs function to estimate flows by birth place, contained the mu element of the returned list.

> # run flow from stock estimation
> library("migest")
> y <- ffs(P1=P1, P2=P2, d=d, b=b)
1 46 
2 0 
> # display with row, col and table totals
> addmargins(y$mu)
, , pob = A

     dest
orig    A   B  C D  Sum
  A   950   0 50 0 1000
  B     0 100  0 0  100
  C     0   0 10 0   10
  D     0   0  0 0    0
  Sum 950 100 60 0 1110

, , pob = B

     dest
orig   A   B  C D Sum
  A   55   0  0 0  55
  B   25 505 25 0 555
  C    0   0 50 0  50
  D    0   0  0 5   5
  Sum 80 505 75 5 665

, , pob = C

     dest
orig   A  B   C  D Sum
  A   80  0   0  0  80
  B   10 30   0  0  40
  C    0  0 800  0 800
  D    0  0   0 40  40
  Sum 90 30 800 40 960

, , pob = D

     dest
orig   A  B C   D Sum
  A   20  0 0   0  20
  B    0 25 0   0  25
  C   10 10 0   0  20
  D   10 10 0 180 200
  Sum 40 45 0 180 265

, , pob = Sum

     dest
orig     A   B   C   D  Sum
  A   1105   0  50   0 1155
  B     35 660  25   0  720
  C     10  10 860   0  880
  D     10  10   0 225  245
  Sum 1160 680 935 225 3000

The fm function returns the flow matrix aggregated over the place of birth dimension in the mu array.

> # display aggregate flows
> f <- fm(y$mu)
> addmargins(f)
     dest
orig   A  B  C D Sum
  A    0  0 50 0  50
  B   35  0 25 0  60
  C   10 10  0 0  20
  D   10 10  0 0  20
  Sum 55 20 75 0 150

….and there you have it, an estimated flow matrix that matches the changes in the stock tables whilst controlling for births and deaths. In the paper I run the code on real migrant stock data provided by the World Bank, to estimate global migrant flow tables.

The ffs function has some different methods to control for deaths in the estimation procedure. The estimation is based on a three way iterative proportional fitting scheme to estimate parameters in a log-linear model, not to dissimilar to that used in a paper based on my Southampton M.Sc. dissertation.

Publication Details:

Abel, G. J. (2013). Estimating global migration flow tables using place of birth data. Demographic Research, 28, 505–546. doi:10.4054/DemRes.2013.28.18

International migration flow data often lack adequate measurements of volume, direction and completeness. These pitfalls limit empirical comparative studies of migration and cross national population projections to use net migration measures or inadequate data. This paper aims to address these issues at a global level, presenting estimates of bilateral flow tables between 191 countries. A methodology to estimate flow tables of migration transitions for the globe is illustrated in two parts. First, a methodology to derive flows from sequential stock tables is developed. Second, the methodology is applied to recently released World Bank migration stock tables between 1960 and 2000 (Özden et al. 2011) to estimate a set of four decadal global migration flow tables. The results of the applied methodology are discussed with reference to comparable estimates of global net migration flows of the United Nations and models for international migration flows. The proposed methodology adds to the limited existing literature on linking migration flows to stocks. The estimated flow tables represent a first-of-a-kind set of comparable global origin destination flow data.

Bank of England Fan Charts in R

I have/will update this post as I expanded the fanplot package.

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 Chart 5.3 in February 2013 Inflation Report

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, in the form of mode, uncertainty and a skewness parameters of a split-normal distribution that underlie their fan charts (The Bank of England predominately refer to the equivalent, re-parametrised, two-piece normal distribution). The probability density of the split-normal distribution is given by Julio (2007) as

f(x; \mu, \sigma_1, \sigma_2) = \left\{  \begin{array}{ll}  \frac{\sqrt 2}{\sqrt\pi (\sigma_1+\sigma_2)} e^{-\frac{1}{2\sigma_1^2}(x-\mu)^2} & \mbox{for } -\infty < x \leq \mu \\  \frac{\sqrt 2}{\sqrt\pi (\sigma_1+\sigma_2)} e^{-\frac{1}{2\sigma_2^2}(x-\mu)^2} & \mbox{for } \mu < x < \infty \\  \end{array},  \right.

where \mu represents the mode parameter, and the two standard deviations \sigma_1 and \sigma_2 can be derived given the overall uncertainty parameter, \sigma and skewness parameters, \gamma, as;

\sigma^2=\sigma^2_1(1+\gamma)=\sigma^2_2(1-\gamma).

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. I added two data objects to the fanplot package to help readers reproduce my plots below. The first, cpi, is a time series object with past values of CPI index. The second, boe, is a data frame with historical details on the split normal parameters for CPI inflation between Q1 2004 to Q4 2013 forecasts by the Bank of England.

> library("fanplot")
> head(boe)
  time0    time mode uncertainty skew
1  2004 2004.00 1.34      0.2249    0
2  2004 2004.25 1.60      0.3149    0
3  2004 2004.50 1.60      0.3824    0
4  2004 2004.75 1.71      0.4274    0
5  2004 2005.00 1.77      0.4499    0
6  2004 2005.25 1.68      0.4761    0

The first column time0 refers to the base year of forecast, the second, time indexes future projections, whilst the remaining three columns provide values for the corresponding projected mode (\mu), uncertainty (\sigma) and skew (\gamma) parameters:

Users can replicate past Bank of England fan charts for a particular period after creating a matrix object that contains values on the split-normal quantile function for a set of user defined probabilities. For example, in the code below, a subset of the Bank of England future parameters of CPI published in Q1 2013 are first selected. Then a vector of probabilities related to the percentiles, that we ultimately would like to plot different shaded fans for, are created. Finally, in a for loop the qsplitnorm function, calculates the values for which the time-specific (i) split-normal distribution will be less than or equal to the probabilities of p.

# select relevant data
y0 <- 2013
boe0 <- subset(boe, time0==y0)
k <- nrow(boe0)

# guess work to set percentiles the BOE are plotting
p <- seq(0.05, 0.95, 0.05)
p <- c(0.01, p, 0.99)
 
# quantiles of split-normal distribution for each probability
# (row) at each future time point (column)
cpival <- matrix(NA, nrow = length(p), ncol = k)
for (i in 1:k) 
   cpival[, i] <- qsplitnorm(p, mode = boe0$mode[i], 
                                sd = boe0$uncertainty[i],
                                skew = boe0$skew[i])

The new object cpival contains the values evaluated from the qsplitnorm function in 6 rows and 13 columns, where rows represent the probabilities used in the calculation p and columns represent successive time periods.

The object cpival can then used to add a fan chart to the active R graphic device. In the code below, the area of the plot is set up when plotting the past CPI data, contained in the time series object cpi. The xlim arguments are set to ensure space on the right hand side of the plotting area for the fan. Following the Bank of England style for plotting fan charts, the background for future values is set to a gray colour, y-axis are plotted on the right hand side, a horizontal line are added for the CPI target and a vertical line for the two-year ahead point.

# past data
plot(cpi, type = "l", col = "tomato", lwd = 2, 
     xlim = c(y0 - 5, y0 + 3), ylim = c(-2, 7), 
     xaxt = "n", yaxt = "n", ylab="")
 
# background
rect(y0 - 0.25, par("usr")[3] - 1, y0 + 3, par("usr")[4] + 1, 
     border = "gray90", col = "gray90")
 
# add fan
fan(data = cpival, data.type = "values", probs = p, 
    start = y0, frequency = 4, 
    anchor = cpi[time(cpi) == y0 - 0.25], 
    fan.col = colorRampPalette(c("tomato", "gray90")),  
    ln = NULL, rlab = NULL)
 
# boe aesthetics
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)
abline(h = 2)  #boe cpi target
abline(v = y0 + 1.75, lty = 2)  #2 year line

boe1
The fan chart itself is outputted from the fan function, where arguments are set to ensure a close resemblance of the R plot to that produced by the Bank of England. The first three arguments in the fan function called in the above code, provide the cpival data to plotted, indicate that the data are a set of calculated values (as opposed to simulations) and provide the probabilities that correspond to each row of cpival object. The next two arguments define the start time and frequency of the data. These operate in a similar fashion to those used when defining time series in R with the ts function. The anchor argument is set to the value of CPI before the start of the fan chart. This allows a join between the value of the Q1 2013 observation and the fan chart. The fan.col argument is set to a colour palette for shades between tomato and gray90. The final two arguments are set to NULL to suppress the plotting of contour lines at the boundary of each shaded fan and their labels, as per the Bank of England style.

Default Fan Chart Plot.

By default, the fan function treats objects passed to the data argument as simulations from sequential distributions, rather than user-created values corresponding probabilities provided in the probs argument (as above). An alternative plot below, based on simulated data and default style settings in the fan function produces a fan chart with a greater array of coloured fans with labels and contour lines alongside selected percentiles of the future distribution. To illustrate we can simulate 10,000 values from the future split-normal distribution parameters from Q1 2013 in the boe0 data frame using the rsplitnorm function

#simulate future values
cpisim <- matrix(NA, nrow = 10000, ncol = k)
for (i in 1:k) 
   cpisim[, i] <- rsplitnorm(n=10000, mode = boe0$mode[i],
                                      sd = boe0$uncertainty[i], 
                                      skew = boe0$skew[i])

The fan chart based on the simulations in cpisim can then be added to the plot;

# truncate cpi series
cpi0 <- ts(cpi[time(cpi)<2013], start=start(cpi), 
           frequency=frequency(cpi) )

# past data
plot(cpi0, type = "l", lwd = 2, 
     xlim = c(y0 - 5, y0 + 3.25), ylim = c(-2, 7))

# add fan
fan(data = cpisim, start = y0, frequency = 4) 

boe2
The fan function calculates the values of 100 equally spaced percentiles of each future distribution when the default data.type = "simulations" is set. This allows 50 fans to be plotted from the heat.colours colour palate, providing a finer level of shading in the representation of future distributions. In addition, lines and labels are provided along each decile. The fan chart does not connect to the last observation as anchor = NULL by default.