\HeaderA{ProbForecastGOP}{Probabilistic Weather Forecasts using the GOP method}{ProbForecastGOP}
\keyword{file}{ProbForecastGOP}
\begin{Description}\relax
This function generates probabilistic forecasts of weather fields 
using the Geostatistical Output Perturbation method as described by Gel, 
Raftery and Gneiting (2004).
\end{Description}
\begin{Usage}
\begin{verbatim}
ProbForecastGOP(day, obs, forecast, id, coord1, coord2, cut.points=NULL, 
max.dist=NULL, nbins=300, variog.model="exponential", max.dist.fit=NULL, 
init.val=NULL, fix.nugget=FALSE, coord1.grid, coord2.grid, forecast.grid, 
n.sim=99, out="SIM", n.displ=4, qt.displ=c(10,50,90))
\end{verbatim}
\end{Usage}
\begin{Arguments}
\begin{ldescription}
\item[\code{day}] numeric vector containing the day of observation.
\item[\code{obs}] numeric vector containing the observed weather quantity.
\item[\code{forecast}] numeric vector containing the forecasted weather quantity.
\item[\code{id}] vector with the id of the metereological stations.
\item[\code{coord1}] vector containing the longitudes of the metereological stations.
\item[\code{coord2}] vector containing the latitudes of the metereological 
stations.
\item[\code{cut.points}] numeric vector containing the cutpoints used for variogram binning.
\item[\code{max.dist}] a numerical value giving the upper bound for the 
distance considered in the variogram computation.
\item[\code{nbins}] a numerical value giving the number of bins for variogram binning.
\item[\code{variog.model}] character string giving the name of the parametric model to fit to the empirical variogram. Implemented models are \kbd{exponential}, 
\kbd{spherical}, \kbd{gauss}, \kbd{gencauchy}, and \kbd{matern}.
\item[\code{max.dist.fit}] number giving the maximum distance considered when fitting the variogram.
\item[\code{init.val}] numeric vector giving the initial values for the parameters in the variogram model. The number of initial values to be entered depends on the parametric model specified. 

If the variog.model specified is \kbd{exponential}, \kbd{spherical} or 
\kbd{gauss}, then the parameters required are, in order, the nugget 
effect, the variance and the range. 

If the variog.model specified is \kbd{gencauchy}, the parameters required 
are, in order, the nugget effect, the variance, the range, the smoothness 
parameter \kbd{a} and the long-range parameter \kbd{b}.

If the variog.model specified is \kbd{matern} the parameters required are, 
in order, the nugget effect, the variance, the range, and the smoothness 
parameter \kbd{a}.

For more details on the valid range for the parameters and for the 
equation 
of the variogram models listed above, look below at the section "Details".
\item[\code{fix.nugget}] logical field indicating whether the nugget should be considered fixed or not. 

If \kbd{TRUE} the nugget effect will be assumed to be constant, and a value for the fixed nugget effect can be also provided. If the value provided is different from the one entered in the \code{init.val} field, then the value of the nugget effect is taken to be the one entered in the \code{init.val} field. 

If \kbd{FALSE} the nugget effect will be estimated along with the other parameters in the variogram model.
\item[\code{coord1.grid}] numeric vector containing the longitudes of the grid points for the forecast field.
\item[\code{coord2.grid}] numeric vector containing the latitudes of the grid points for the forecast field.
\item[\code{forecast.grid}] numeric vector containing the forecast grid.
\item[\code{n.sim}] number of realizations to be simulated.
\item[\code{out}] character string indicating the output produced by the function. Possible values are \kbd{VARIOG}, \kbd{FIT} and \kbd{SIM}. 

If \kbd{VARIOG} is specified, the function computes the empirical 
variogram of the forecast errors averaged over time and displays a plot of the empirical variogram. 

If \kbd{FIT} is specified, the function, in addition, fits a parametric model to the empirical variogram, returns estimates of the parameters and display both the empirical and the fitted variogram.

If \kbd{SIM} is specified, the function simulates also realizations of the weather random field using the parametric variogram model specified. The function displays the simulated weather fields and some specified percentiles.
\item[\code{n.displ}] number of realizations to be displayed on screen.
\item[\code{qt.displ}] numeric vector containing the quantiles to be displayed.
\end{ldescription}
\end{Arguments}
\begin{Details}\relax
This function generates probabilistic weather field forecasts by generating ensemble members using the Geostatistical Output Perturbation method presented by Gel, Raftery and Gneiting (2004).

The Geostatistical Output Perturbation method consists of the following steps:

1. determination of the residuals, by regressing the forecasts on the observed weather quantity and computation of their empirical variogram. 

2. estimation of the parameters of the variogram model fitted to the empirical variogram. 

3. simulation of weather random fields using the estimated parametric variogram model.


- \emph{Computation of the empirical variogram of the residuals}

The empirical variogram of the residuals is calculated by determining, for each day, the distance among all pairs of stations that have been observed in the same day and by calculating for each day the sum of all the squared differences in the residuals within each bin. These sums are then averaged over time, with weights for each bin given by the sum over time of the number of pairs of stations within the bin. 

The formula used is:
\deqn{ \gamma(h) = \sum_d \frac{1}{2N_{(h,d)}} (\sum_h (Y(i+h,d)-Y(i,d))^2)}{}


- \emph{Estimation of the parameters of the parametric variogram model 
fitted to the empirical variogram.}

The function estimates the parameters in the variogram model by minimizing the weighted least-square loss function

\deqn{LOSS(\theta) = \sum_k n_k [ \left( \frac{\hat{\gamma_k} - \gamma_k(\theta)}{\gamma_k(\theta)} \right)^2 ]}{}

where \eqn{n_k}{} is the number of pairs contributing to the variogram 
computation in the \eqn{k}{}th bin, \eqn{\hat{\gamma_k}}{} is the value of the 
empirical variogram in the \eqn{k}{}th bin and \eqn{\gamma_k(\theta)}{} is the value of the estimated parametric variogram model at the midpoint of the \eqn{k}{}th bin.

- Parametric variogram models -

The parametric variogram models implemented are: exponential, spherical, 
gaussian, generalized Cauchy and Whittle-Matern.

- \kbd{exponential}

The equation of the exponential variogram with parameters the nugget 
effect, \eqn{\rho}{}, the variance, \eqn{\sigma^{2}}{}, and the range, 
\eqn{r}{}, is given by:

\deqn{\gamma(d) = \rho+\sigma^{2} \cdot (1-exp(- \frac{d}{r}))}{}

where \eqn{d}{} is the distance between the two locations and 
\eqn{\gamma(d)}{} is the value of the exponential variogram at distance 
\eqn{d}{}. The nugget effect, \eqn{\rho}{}, is a non-negative number, while 
the variance, \eqn{\sigma^2}{}, and the range, \eqn{r}{}, are positive 
numbers.


- \kbd{spherical}

The equation of the spherical variogram with parameters the nugget effect, 
\eqn{\rho}{}, the variance, \eqn{\sigma^2}{}, and the range, \eqn{r}{}, is given 
by:

\deqn{\gamma(d) = \rho+\sigma^{2} \cdot (\frac{3}{2} 
\cdot\frac{d}{r}-\frac{1}{2} \cdot \frac{d^3}{r^3})}{}

where \eqn{d}{} is the distance between the two locations and 
\eqn{\gamma(d)}{} is the value of the spherical variogram at distance 
\eqn{d}{}. The nugget effect, \eqn{\rho}{}, is a non-negative number, while 
the variance, \eqn{\sigma^2}{}, and the range, \eqn{r}{}, are positive 
numbers.


- \kbd{gauss}

The equation of the gaussian variogram with parameters the nugget effect, 
\eqn{\rho}{}, the variance, \eqn{\sigma^2}{}, and the range, \eqn{r}{}, is given 
by:

\deqn{\gamma(d) = \rho+\sigma^{2} \cdot (1-exp(- \frac{d^2}{r^2} ))}{}

where \eqn{d}{} is the distance between the two locations and 
\eqn{\gamma(d)}{} is the value of the gaussian variogram at distance 
\eqn{d}{}. The nugget effect, \eqn{\rho}{}, is a non-negative number, while 
the variance, \eqn{\sigma^2}{}, and the range, \eqn{r}{}, are positive 
numbers.


- \kbd{gencauchy}

The equation of the generalized Cauchy variogram with parameters the 
nugget effect, \eqn{\rho}{}, the variance, \eqn{\sigma^2}{}, the range, 
\eqn{r}{}, 
the smoothness parameter, \eqn{a}{}, and the long-range parameter, \eqn{b}{}, is given by:

\deqn{\gamma(d) = \rho+\sigma^{2} \cdot 
(1-(1+\frac{d^a}{r^a})^{- \frac{b}{a}})}{}

where \eqn{d}{} is the distance between the two locations and 
\eqn{\gamma(d)}{} is the value of the generalized Cauchy variogram at 
distance \eqn{d}{}. The nugget effect, \eqn{\rho}{} is a non-negative number, 
the variance, \eqn{\sigma^2}{} and the range, \eqn{r}{}, are positive numbers, 
the smoothness parameter, \eqn{a}{}, is a number in \eqn{(0,2]}{}, and the long-range parameter, \eqn{b}{}, is a positive number.


- \kbd{matern}

The equation of the Whittle-Matern variogram with parameters the nugget 
effect, \eqn{\rho}{}, the variance, \eqn{\sigma^2}{}, the range, \eqn{r}{}, and 
the smoothness parameter, \eqn{a}{}, is given by:

\deqn{\gamma(d) = \rho+\sigma^{2} \cdot (1-(\frac{2^{1-a}}{\Gamma(a)}\cdot 
\frac{d^a}{r^a} \cdot K_{a}(\frac{d}{r}))}{}

where \eqn{d}{} is the distance between the two locations, \eqn{\gamma(d)}{} 
is the value of the Whittle-Matern variogram at distance \eqn{d}{}, 
\eqn{\Gamma}{} is the gamma function and \eqn{K_a}{} is the Bessel function of 
the third kind with parameter \eqn{a}{}. The nugget effect, \eqn{\rho}{}, is a 
non-negative number, the variance, \eqn{\sigma^2}{}, the range, \eqn{r}{}, and the smoothness parameter, \eqn{a}{}, are positive numbers.


- Defaults-

If the vector with the cutpoints is not specified, the cutpoints are determined so that there are \code{nbins} bins with approximately the same number of pairs per bin. If both the vector with the cutpoints and the number of bins, \code{nbins}, are unspecified, the function by default determines the cutpoints so that there are 300 bins with approximately the same number of pairs per bin. 

The default value for the maximum distance considered in the variogram computation is the \eqn{90}{}-th percentile of the distances between all stations. 

The default value for the maximum distance is \eqn{\frac{1}{2 \cdot 
\sqrt{2}}}{} times the maximum distance considered in the empirical 
variogram. 

Default for the initial values of the parameters is \kbd{NULL}, in which 
case the initial values for the parameters are determined internally and depend on the empirical variogram values. By default, the nugget effect is not considered constant. Thus, the default value for the \code{fix.nugget} field is \kbd{FALSE}.

The default for \code{out} is \kbd{SIM}: by default, \eqn{99}{} weather 
random fields are simulated, and \eqn{4}{} of them are displayed. If no 
vector of percentiles, is provided, the \eqn{10}{}th, \eqn{50}{}th and \eqn{90}{}th percentiles are determined and displayed.
\end{Details}
\begin{Value}
The function returns both a numerical and a graphical summary. The graphical and numerical summary provided depend on the specification given in the field \code{out}.

The graphical summary consists of a plot of the empirical variogram, if 
\code{out} is equal to \kbd{VARIOG}. 

It consists of a plot of the empirical variogram with superimposed the fitted parametric model, if \code{out} is equal to \kbd{FIT}. 

If \code{out} is equal to \kbd{SIM}, the graphical summary consists of 
several plots: a plot of the empirical variogram with superimposed the fitted 
parametric model, plots of the simulated weather fields, and plots of 
specified percentiles of the weather random field. The plots are printed on multiple pages, and before displaying each page, the user will be asked for input.

The numerical summary include the following:
\begin{ldescription}
\item[\code{bias.coeff}] additive and multiplicative bias obtained by linear regression of the forecasted weather quantity on the observed weather quantity. Always provided.
\item[\code{se.bias.coeff}] standard errors of the bias coefficients. Always provided.
\item[\code{res.var}] variance of the forecast errors. Always provided.
\item[\code{bin.midpoints}] numeric vector with midpoints of the bins used in the empirical variogram computation. Always provided.
\item[\code{number.pairs}] numeric vector with the number of pairs per bin. Always provided.
\item[\code{empir.variog}] numeric vector with the empirical variogram values. Always provided.
\item[\code{model}] name of the parametric model fitted to the empirical variogram. Provided only if \code{out} is equal to \kbd{FIT} or \kbd{SIM}.
\item[\code{nugget}] Estimate of the nugget effect. Provided only if \code{out} is equal to \kbd{FIT} or \kbd{SIM}.
\item[\code{variance}] Estimate of the variance. Provided only if \code{out} is equal to \kbd{FIT} or \kbd{SIM}.
\item[\code{range}] Estimate of the range. Provided only if \code{out} is 
equal to \kbd{FIT} or \kbd{SIM}.
\item[\code{additional.par}] Numeric vector with the estimates of the 
additional parameters required by the parametric variogram model (only if 
the model specified is either \kbd{whittlematern} or \kbd{gencauchy}). 
This is provided only if \code{out} is equal to \kbd{FIT} or \kbd{SIM}.
\item[\code{sim.fields}] 3-dimensional array where each layer contains the values of the simulated weather field at the gridded locations. Provided only if \code{out} is equal to \kbd{SIM}.
\item[\code{pct.fields}] 3-dimensional array where each layer contains the specified percentiles of the weather field. Provided only if \code{out} is equal to \kbd{SIM}.
\end{ldescription}
\end{Value}
\begin{Note}\relax
Computation of the empirical variogram might require substantial 
computing time. Hence, it is advisable to run the function once with 
\code{out} equal to \kbd{VARIOG}, save part of the output and use it as input for the \code{Variog.fit} function to estimate a parametric variogram model. Finally, \code{Field.sim} will simulate realizations of the weather field from the specified geostatistical model.
\end{Note}
\begin{Author}\relax
Gel, Y., Raftery, A. E., Gneiting, T., Berrocal, V. J. \email{veronica@stat.washington.edu}.
\end{Author}
\begin{References}\relax
Gel, Y., Raftery, A. E., Gneiting, T. (2004). Calibrated 
probabilistic mesoscale weather field forecasting: The Geostatistical 
Output Perturbation (GOP) method (with discussion). \emph{Journal of the 
American Statistical Association}, \bold{Vol. 99 (467)}, 575--583.

Gel, Y., Raftery, A. E., Gneiting, T., Berrocal, V. J. (2004). 
Rejoinder. \emph{Journal of the American Statistical Association}, 
\bold{Vol. 99 (467)}, 588--590.  

Cressie, N. A. C. (1993). \emph{Statistics for Spatial Data} (revised ed.). Wiley: New York.

Gneiting, T., Schlather, M. (2004). Stochastic models that separate 
fractal dimension and the Hurst effect. \emph{SIAM Review}, \bold{46}, 269--282.

Stein, M. L. (1999). \emph{Interpolation of Spatial Data - Some Theory for Kriging}. Springer-Verlag: New York.

Schlather, M. (2001). Simulation and analysis of Random Fields. \emph{R 
News 1(2)}, 18--20.

Nychka, D. (2004). \emph{The fields package}. Available at: http:lib.stat.cmu.edu/R/CRAN/doc/package/fields.pdf.
\end{References}
\begin{SeeAlso}\relax
\code{\LinkA{Emp.variog}{Emp.variog}} for computation of empirical variogram 
of forecast errors averaged over time, \code{\LinkA{Variog.fit}{Variog.fit}} for 
estimation of parameters in a geostatistical model using Weighted Least 
Square, \code{\LinkA{Field.sim}{Field.sim}} for realizations of weather fields with 
the specified covariance structure.
\end{SeeAlso}
\begin{Examples}
\begin{ExampleCode}
## Loading data
library(fields)
library(RandomFields)
data(slp)
day <-slp$date.obs
id <- slp$id.stat
coord1 <- slp$lon.stat
coord2 <- slp$lat.stat
obs <- slp$obs
forecast <-slp$forecast

data(gridlong)
coord1.grid <- gridlong$gridded.lon

data(gridlat)
coord2.grid <- gridlat$gridded.lat

data(forecast.grid)
forecast.grid <- forecast.grid$gridded.forecast

## Specified cutpoints, default values for all the other fields.
## Only empirical variogram computation
empirical.variog <- ProbForecastGOP(day=day,obs=obs,forecast=forecast,
id=id,coord1=coord1,coord2=coord2,cut.points=seq(0,1000,by=2),
max.dist=NULL,nbins=NULL,variog.model="exponential",max.dist.fit=NULL,
init.val=NULL,fix.nugget=FALSE,coord1.grid=coord1.grid,
coord2.grid=coord2.grid,forecast.grid=forecast.grid,n.sim=10,
out="VARIOG",n.displ=4,qt.displ=c(5,50,95))

## Unspecified cutpoints.
## Fit of a parametric variogram to an empirical variogram.
fitted.variog <- ProbForecastGOP(day=day,obs=obs,forecast=forecast,id=id,
coord1=coord1,coord2=coord2,cut.points=NULL,max.dist=NULL,nbins=NULL,
variog.model="exponential",max.dist.fit=NULL,init.val=NULL,fix.nugget=FALSE,
coord1.grid=coord1.grid,coord2.grid=coord2.grid,forecast.grid=forecast.grid,
n.sim=10,out="FIT",n.displ=4,qt.displ=c(5,50,95))

## Unspecified cutpoints. 
## Whole routine.
simulation.fields <- 
ProbForecastGOP(day=day,obs=obs,forecast=forecast,id=id,coord1=coord1,
coord2=coord2,cut.points=NULL,max.dist=NULL,nbins=NULL,variog.model=NULL,
max.dist.fit=NULL,init.val=NULL,fix.nugget=FALSE,coord1.grid=coord1.grid,
coord2.grid=coord2.grid,forecast.grid=forecast.grid,n.sim=4,out="SIM",
n.displ=4,qt.displ=c(5,50,95))


\end{ExampleCode}
\end{Examples}

