This file explains how to use the R code in "hierarchical.txt" to fit the model discussed in "A Hierarchical Model for Extreme Wind Speeds"  to the wind speed data collected by the U.K. Meteorological Office in northern/central England.  

The file "hierarchical.txt" should be copy-and-pasted into an R terminal.  
The 12 sets of threshold exceedances should also be in the current working directory. 

 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Section 1 of the code reads the data into R, and separates the data from each site into its monthly (seasonal) components (e.g. "janshef.exceed" is the set of threshold exceedances in January at Sheffield; "aug4238.exceed" is the set of threshold exceedances in August at Finningley). 
     
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Section 2 of the code contains the function which is used to fit the hierarchical model.  The function is called "random.effects" and has 13 arguments:

n:	      the number of iterations in the MCMC scheme
seasons:      a combined vector which contains all threshold excesses by season - 
	      e.g. "all.months" 
sites:	      a combined vector which contains all threshold excesses by site -
	      e.g. "all.sites"
tau.etastart: a starting value for the parameter \tau_{\tilde{\sigma}}
tau.xistart:  a starting value for the parameter \tau_{\xi}
zeta.etastart:a starting value for the parameter \zeta_{\tilde{\sigma}}
zeta.xistart: a starting value for the parameter \zeta_{\xi}
bstart-gstart:hyper-parameters for the means and variances of the random effects 
	      parameters

(Throughout, the parameter "eta" is used in place of "log(\tilde{\sigma})")

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In Section 2.1 we specify the starting values for the random effects parameters themselves - e.g. "gamma.eta1start" is the starting value for \gamma_{\tilde{\sigma}}^{(1)}, i.e. the seasonal effect for January for the GPD scale parameter.  

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In Section 2.2 we specify the variances of the random walk innovations for each of the random effects parameters.   

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In Section 2.3 we set up the vectors which will take the candidate values at each step in the MCMC scheme (e.g. "cangamma.eta1" will take the candidate values for the parameter \gamma_{\tilde{\sigma}}^{(1)}).  We also specify the vectors which will take the posterior draws for each of the random effects parameters.  We use:

x: for the seasonal effects for the GPD scale parameter (\gamma_{\tilde{\sigma}})
y: for the seasonal effects for the GPD shape parameter (\gamma_{\xi})
v: for the site effects for the GPD scale parameter (\epsilon_{\tilde{\xi}})
w: for the site effects for the GPD shape parameter (\epsilon_{\xi})

The numbers 1, ..., 12 after each of the x, y, w, v correspond to the season or site to which the posterior draws belong (for seasons, 1=Jan, ..., 12=Dec; for sites, 1=Sheffield, ..., 12= Ringway (Manchester))

Vectors such as "eta.one.one" are the recovered posterior draws for each of the GPD parameters (in this case log(\tilde{\sigma}) for season 1, site 1), found by adding together the corresponding seasonal and site effects.  

Vectors labelled "aprob*" will contain the acceptance probabilities for each of the random effects parameters in the MCMC scheme.  

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In Section 2.4 we specify the GPD log-likelihood function in terms of the seasonal and site effects for the GPD parameters and the threshold vector U.  

The vector "thresh" contains the threshold used to identify extreme wind speeds.  This vector is of length 144, and consists of the twelve monthly thresholds for wind speeds at Sheffield, Bradfield, ..., Ringway (Manchester).  These have been obtained from mean residual life plots.  

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Section 2.5 starts the MCMC iterations.  Each random effects parameter has been initialised at a starting value, and so the iterations run from i=1 up to n, where n is the total number of iterations.  A Metropolis-within-Gibbs scheme is used, sampling from the full conditionals where these are specified, and using a Metropolis random walk update elsewhere.  Each of the parameters \gamma_{\tilde{\sigma}}^{(1)}, \gamma_{xi}^{(1)}, \epsilon_{\tilde{\sigma}}^{(1)}, \epsilon_{\xi}^{(1)}, ...,\gamma_{\tilde{\sigma}}^{(12)}, \gamma_{xi}^{(12)}, \epsilon_{\tilde{\sigma}}^{(12)}, \epsilon_{\xi}^{(12)} are updated in the MCMC scheme at each iteration.  

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In Section 2.6 we specify what the program should return - all acceptance probabilities and draws from the posteriors for each of the random effects parameters.  


NOTE 1
The file "hierarchical.txt" does not contain the additional code for obtaining posterior draws for the site effects for the logistic dependence parameter \alpha.  This is done separately using the posterior draws for the GPD parameters to transform the marginal distributions for each site to unit Frechet, and then implementing a Metropolis-Hastings algorithm to obtain corresponding posterior draws for \alpha using the Markov chain model likelihood in equation (2) to compute acceptance probabilities.  For more information on how this is done, see Fawcett (2005). 

NOTE 2
For more help on how to implement this code in R, contact Lee Fawcett: lee.fawcett@ncl.ac.uk