Demo file for the fanplot package

I have added a demo file to the latest version of the fanplot package. It has lots of examples of different plotting styles to represent uncertainty in time series data. In the updated package I have added functionality to plot fan charts based on irregular time series objects from the zoo package, plus the use of alternative colour palettes from the RColorBrewer and colorspace packages. All plots are based on the th.mcmc object, the estimated posterior distributions of the volatility in daily returns from the Pound/Dollar exchange rate from 02/10/1981 to 28/6/1985. To run the demo file from your R console (ensure fanplot, zoo, tsbugs, RColorBrewer and colorspace packages are all installed beforehand);

# if you want plots in separate graphic devices 
# do not run this first line...
par(mfrow = c(10,2))
# run demo
demo(sv_fan, package = "fanplot", ask = FALSE)

The demo script should output this set of plots:
If you wish, click on the image above and take a closer look in your browser. In R, you can save the PDF version of all the plots on one graphics device (which looks much better than what comes up in my R graphics device):

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

You can also view the demo file for a closer look at the arguments used in each plot:"demo/sv_fan.R", package = "fanplot"))

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;


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

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) 

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.

The tsbugs package for R

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

My tsbugs package has gone up on CRAN. I decided not to write a vignette for the submission, as it would have involved doing some estimation in BUGS via R2WinBUGS or R2OpenBUGS and running into some problems when submitted the package. Instead I thought I would post a quick guide below. Note, the package code is on my github if you would like to contribute more models to the package…

The functions in the tsbugs package are aimed to automate the writing of time series models to run in WinBUGS or OpenBUGS. I created these functions a while back when I was doing some work on model averaging for time series models. I found it a lot easier to build R functions to write the BUGS models than the more error-inducing process of copy and pasting BUGS scripts, and then making slight alterations to create new models. It also allowed me to add arguments to specify different lag lengths, prior distributions, variance assumptions and data lengths. Below are examples for three types of time series models; autorgressive models with constant variance, stochastic volatility and random variance shift models.

Autoregressive Models

The ar.bugs command builds a BUGS script for autoregressive (AR) models ready to use in R2OpenBUGS. For example, consider the LakeHuron data.

LH <- LakeHuron 
plot(LH, main="Level (ft)") 
plot(diff(LH), main="Differenced Level")

We can construct a AR(1) model for this data (after differencing the data to obtain a stationary mean) as such:

ar1 <- ar.bugs(y=diff(LH), ar.order=1)

The ar.bugs function allows for alternative specifications for prior distributions, forecasts and the inclusion of mean term:

ar2 <- ar.bugs(y=diff(LH), ar.order=2, ar.prior="dunif(-1,1)", 
               var.prior="dgamma(0.001,0.001)", k = 10, 
               mean.centre = TRUE)

The tsbugs objects can be used with R2OpenBUGS to easily run models from R. This is made even easier using the inits and nodes functions (also in the tsbugs package). For example:

writeLines(ar2$bug, "ar2.txt") 
ar2.bug <- bugs(data = ar2$data,
                inits = list(inits(ar2)),
                param = c(nodes(ar2, "prior")$name, ""),
                model = "ar2.txt", 
                n.iter = 11000, n.burnin = 1000, n.chains = 1)

Note, 1) the model is written to a .txt file (as required by R2OpenBUGS), 2) the data used is part of the tsbugs object. The ar.bugs command cleans the data and adds missing values at the end of the series for foretasted values, 3) the initial values offered by the inits function are very crude, and with more complicated data or models, users might be better off specifying there own list of initial values. The parameter traces and posterior distributions can be plotted using the coda package:

param.mcmc <- as.mcmc(ar2.bug$sims.matrix[,nodes(ar2, "prior")$name])

The fanplot package can be used to plot the entire series of posterior predictive distributions. We may also plot (after deriving using the diffinv function) the posterior predictive distributions of the lake level:

# derive future level
ynew.mcmc <- ar2.bug$sims.list$
lhnew.mcmc <- apply(ynew.mcmc, 1, diffinv, xi = tail(LH,1))
lhnew.mcmc <- t(lhnew.mcmc[-1,])

# plot differenced
plot(diff(LH) ,xlim=k0+c(-50,10), main="Differenced Level")

# add fan
k0 <- end(LH)[1]
fan(ynew.mcmc, start=k0+1, rcex=0.5)

# plot undifferenced
plot(LH, xlim=k0+c(-50,10), main="Level")
fan(lhnew.mcmc, start=k0+1, rcex=0.5)


Stochastic Volatility Models

The sv.bugs command builds a BUGS script for stochastic volatility SV models ready to use in R2OpenBUGS. For example, consider the svpdx data.

# plot
plot(svpdx$pdx, type = "l", 
     main = "Return of Pound-Dollar exchange rate data 
             from 2nd October 1981 to 28th June 1985", 
     cex.main = 0.8)

We can construct a AR(0)-SV model for this data, and also obtain posterior simulations using the sv.bugs command:

y <- svpdx$pdx
sv0 <- sv.bugs(y, sim=TRUE)

This model closely matches those presented in
Meyer and Yu (2002). There are further options in the tsbugs package to incorporate different priors that do not involve transformations such as those for psi1 above. Using R2OpenBUGS we can fit the model,

# decent initial value for variance in first period
init <- inits(sv0, warn=FALSE)
init$psi0 <- log(var(y))
# write bug
writeLines(sv0$bug, "sv0.txt")
# might take a while to compile
sv0.bug <- bugs(data = sv0$data,
                inits = list(init),
                param = c(nodes(sv0, "prior")$name,"y.sim","h"),
                model = "sv0.txt", 
                n.iter = 11000, n.burnin = 1000, n.chains = 1)

The volatility and estimates can be easily extracted,

h.mcmc <- sv0.bug$sims.list$h

Which allows us to directly view the estimated volatility process or the time-dependent standard deviation using the fanplot package,

# plot
plot(NULL, xlim = c(1, 945)+c(0,40), ylim = c(-4,2), 
     main="Estimated Volatility from SV Model") 
# fan
fan(h.mcmc, type = "interval")

We can also plot the posterior simulations from the model:

# derive percentiles
y.mcmc <- sv0.bug$sims.list$y.sim
# plot
plot(NULL, type = "l", xlim = c(1, 945)+c(0,20), 
     ylim = range(y), 
     main = "Posterior Model Simulations and Data")


Random Variance Shift Models

The rv.bugs command builds a BUGS script for random variance (RV) shift models, similar to that of McCulloch and Tsay (1993) ready to use in R2OpenBUGS. Consider the ew data.

r <- ts(ew[2:167]/ew[1:166]-1, start=1841)
y <- diff(r) 
plot(y, main="Difference in England and Wales 
              Population Growth Rate")

We can create a BUGS script to fit a RV model to this data, including posterior simulations, using the rv.bugs command:

rv0<-rv.bugs(y, sim=TRUE)

and then run the script in R2OpenBUGS (this can take a couple of hours):

# decent inital value for variance in first period
init <- inits(rv0, warn=FALSE)
# write bug
# might take a while to compile
rv0.bug <- bugs(data = rv0$data,
                inits = list(init),
                param = c(nodes(rv0, "prior")$name,"y.sim",
                model = "rv0.txt", 
                n.iter = 11000, n.burnin = 1000, n.chains = 1)

We can plot the posterior simulations from the model using the fanplot package:

# derive percentiles
y0 <- tsp(y)[1]
y.mcmc <- rv0.bug$sims.list$y.sim

# plot
plot(NULL, xlim=tsp(y)[1:2]+c(-5,5), ylim = range(y), 
     main="Posterior Simulations")
fan(y.mcmc, start = y0, rlab=c(10,50,90), llab=TRUE)

Alongside the posterior distributions of the standard deviations,

# derive sigma
h.mcmc <- rv0.bug$sims.list$h
sigma.mcmc <- sqrt(exp(h.mcmc))
# plots
plot(NULL, xlim =tsp(y)[1:2]+c(-5,5), ylim = c(0,0.008), 
     main="Standard Deviation")
fan(sigma.mcmc, start = y0, rlab=c(5,50,95), llab = c(5,50,95))

The posterior distributions of the probability of a variance shift and multiplier effect of the shift in variance (delta[t] and beta[t] in the BUGS model) can also be plotted. Note, when there is no variance shift, the posterior of the beta[t] is similar to the prior distribution.

#extract data
delta.mcmc <- rv0.bug$sims.list$delta
beta.mcmc <- rv0.bug$sims.list$beta

# plots
plot(NULL, xlim = tsp(y)[1:2]+c(-5,5), ylim = c(0,1),
     main="Probability of Variance Change Point")
fan(delta.mcmc, start=y0, ln = NULL, rlab = NULL)

plot(NULL, xlim = tsp(y)[1:2]+c(-5,5), ylim = c(-2,2), 
     main="Variance Multiplier")
fan(beta.mcmc, start=y0)


The fanplot package for R

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

My fanplot package has gone up on CRAN. Below is a online version of the package vignette…

Visualising Time Series Model Results

The fanplot package can also be used to display uncertainty in estimates from time series models. To illustrate, the packages th.mcmc data frame object contains posterior density distributions of the estimated volatility of daily returns y_t from the Pound/Dollar exchange rate from 02/10/1981 to 28/6/1985. These distributions are from a MCMC simulation from a stochastic volatility model given in Meyer and Yu (2002) where it assumed;

y_t | \theta_t = \exp\left(\frac{1}{2}\theta_t\right)u_t \qquad u_t \sim N(0, 1) \qquad t=1,\ldots,n.

The latent volatilities \theta_t, which are unknown states in a state-space model terminology, are assumed to follow a Markovian transition over time given by the state equations:

\theta_t | \theta_{t-1}, \mu, \phi, \tau^2 = \mu + \phi \log \sigma^2_{t-1} + v_t \qquad v_t \sim N(0, \tau^2) \qquad t=1,\ldots,n

with \theta_0 \sim N(\mu, \tau^2).

The th.mcmc object consists of (1000) rows corresponding to MCMC simulations and (945) columns corresponding to each t. A fan chart of the evolution of the distribution of \theta_t can be visualised using the fanplot package via,

# empty plot
plot(NULL, main="Percentiles", xlim = c(1, 965), ylim = c(-2.5, 1.5))

# add fan
fan(data = th.mcmc)

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 palette, providing a fine level of shading. In addition, lines and labels are provided along each decile.

Prediction Intervals

When argument type = "interval" is set, the probs argument corresponds to prediction intervals. Consequently, the fan chart comprises of 3 different shades, running from the darkest shade for the 50th prediction interval to the lightest for the 95th prediction interval.

# empty plot
plot(NULL, main="Prediction Intervals",
     xlim = c(-20, 965), ylim = c(-2.5, 1.5))

# add fan
fan(data = th.mcmc, type = "interval", llab=TRUE, rcex=0.6)

Contour lines are overlayed for the upper and lower bounds of each prediction intervals, as set using the ln command. A further line is plotted along the median of \theta_t, which is controlled by the med.ln argument (set to TRUE by default when type="interval"). The default labels on the right hand side correspond to the upper and lower bounds of each plotted line. The left labels are added by setting llab = TRUE. Note, some extra room is created for the labels by setting the xlim = c(-20, 965) argument of plotting area to a wider range than the original data (945 observations). The text size of the right labels are controlled using the rcex argument. The left labels, by default, take the same text size as rcex although they can be separately controlled using the lcex argument.

Alternative Colours

Alternative colour schemes to the default heat.colors, can be obtained by supplying a colorRampPalette to the fan.col argument. For example, a new palette running from blue to white, via grey can be passed using;

# empty plot
plot(NULL, main="Alternative Colour Scheme",
xlim = c(-20, 965), ylim = c(-2.5, 1.5))

# add fan
fan(data = th.mcmc, rlab=seq(20,80,15), llab=c(10,50,90),
    fan.col=colorRampPalette(c("royalblue", "grey", "white")))

Alternative labels are specified using the rlab and llab arguments.

Spaghetti Plots

Spaghetti plots can be used to represent uncertainty shown by a range of possible future trajectories or past estimates. For example using the th.mcmc object, 20 random sets of \theta_t can be plotted when setting the argument style = "spaghetti";

# empty plot
plot(NULL, main="Spaghetti Plot", xlim = c(-20, 965), ylim = range(th.mcmc))

# transparent fan with visible lines
fan(th.mcmc, ln=c(5, 50, 95), llab=TRUE, alpha=0, ln.col="orange")

# spaghetti lines
fan(th.mcmc, style="spaghetti", n.spag=20)

The spaghetti lines are superimposed on a fan chart in order to illustrate some underlying probabilities. The initial fan chart is completely transparent from setting the transparency argument alpha to 0. In order for the percentile lines to be visible a non-transparent colour is used for the ln.col argument.

Forecast Fans

The fanplot package can also be used to illustrate probabilistic forecasts. For example, using the auto.arima function in the forecast package a model for the time series for net migration to the United Kingdom (contained in the ips data frame of the fanplot package) can be fitted.

#create time series
net <- ts(ips$net, start=1975)
#fit model
m <- auto.arima(net)
Series: net
ARIMA(1,1,2) with drift         

          ar1      ma1      ma2   drift
      -0.2301  -0.0851  -0.6734  6.7625
s.e.   0.3715   0.3620   0.1924  1.4154

sigma^2 estimated as 1231:  log likelihood=-179.3
AIC=368.6   AICc=370.54   BIC=376.66

We may then simulate 1000 values from the selected model using the simulate.Arima function, and plot the results.

mm <- matrix(NA, nrow=1000, ncol=5)
for(i in 1:1000)
  mm[i,] &amp;lt;- simulate(m, nsim=5)

# empty plot
plot(net, main="UK Net Migration", xlim=c(1975,2020), ylim=c(-100,300))

# add fan
fan(mm, start=2013)

Users might want to connect the fan with the past data. This can be achieved by providing the last value to the anchor argument.

# empty plot
plot(net, main="UK Net Migration",
     xlim=c(1975,2020), ylim=c(-100,300))

# add fan
fan(mm, start=2013, anchor=net[time(net)==2012],
    type="interval", probs=seq(5, 95, 5), ln=c(50, 80))

More shades for the fan are added to the plot (over the default 3 used for a interval fans) by supplying a sequence to the probs argument. Alternative contour lines (from the default median, 50th, 80th and 95th percentiles for interval fans) are added using the ln argument.