fa<-function(x,y,reduction=F,nmodes=3){

# Bayesian forecast assimilation for cross-validated calibration and combination
# of forecasts
#
# Description: This function calls functions faprior, falikelihood, 
# and faprediction and estimates the mean (ya) and the
# square root of the diagonal elements of the
# covariance matrix (d) of the posterior distribution
# y|x ~ N(ya,d)
#
# Usage: fa(x,y,reduction,nmodes)
#
# Arguments:
# x:    n x p data matrix containing time series of p forecast variables
# y:    n x q data matrix containing time series of q observable
#	variables
# reduction: Logical. If TRUE performs maximum covariance analysis
#            on the cross-covariance matrix (y'*x) and uses the
#            selected number of modes ('nmodes', see below) to model the
#            likelihood function. Default is FALSE.
# nmodes: number of modes used in the likelihood modelling. Default
#         is 3. Argument only used if 'reduction' is TRUE.  
#
# Outputs:
# $ya: q x n matrix containing q predicted observable state time series
# $d : q x q matrix containing predicted observable error variance
#      (i.e. the square root of the diagonal elementes of observable 
#       error covariance matrix d.)
#
# Authors: Caio Coelho <c.a.d.s.coelho@reading.ac.uk>       13 Feb 2006
#          David Stephenson <d.b.stephenson@reading.ac.uk>
#
# Reference: Stephenson D. B., C.A.S. Coelho, F. J. Doblas-Reyes and
#            M. Balmaseda, 2005: Forecast Assimilation: A unified 
#            framework for the combination of multi-model weather and
#            climate predictions.  Tellus A . Vol.  57, 253-264.
#
# Example 1: Calibration and combination of forecasts for a single location
#
#          xfa<-read.table("xfa.txt") # xfa.txt contains 3 time series of forecasts 
#          yfa<-read.table("yfa.txt") # yfa.txt contains 1 time series of observations
#          out<-fa(xfa,yfa) 
#          out$ya # combined and calibrated forecast mean 
#          out$d  # combined and calibrated forecast standard deviation
#
# Example 2: Calibration and combination of grided coupled model
#            forecasts
#
#          xfa<-read.table("xfasst.txt") # xfa.txt contains 168 time series of forecasts  
#          # Columns 1-56 of xfa.txt are forecasts produced by model A
#          # Columns 57-112 of xfa.txt are forecasts produced by model B
#          # Columns 113-168 of xfa.txt are forecasts produced by model C
#          yfa<-read.table("yfasst.txt") # yfa.txt contains 56 time series of observations
#          out<-fa(xfa,yfa,reduction=T,3) 
#          out$ya # combined and calibrated forecast mean
#          out$d  # combined and calibrated forecast standard deviation


y <- as.matrix(y)
x <- as.matrix(x)
n <- dim(x)[1]
ya<-matrix(nrow=n,ncol=ncol(y))
d <-matrix(nrow=n,ncol=ncol(y))
if(!reduction){
# Cross-validation loop
nn <- 1
while(nn < n + 1) {
xcross <- as.matrix(x[-nn,])	
ycross <- as.matrix(y[-nn,])
# Estimate prior parameters yb (1 x q) and c (q x q)
prior <- faprior(ycross)
# Estimate likelihood parameters g (p x q), yzerogt (n-1 x p) and
# s (p x q)
lik <- falikelihood(xcross, ycross)
# Estimate posterior parameters ya (1 x q) and d (q x q)
out <- faprediction(x[nn,], ycross, prior$yb, prior$c, lik$g, lik$yzerogt,
lik$s)
ya[nn,]<-out$ya
d[nn,]<-sqrt(diag(out$d))
nn <- nn + 1
} # end of cross-validation loop
}
else{
# Cross-validation loop
nn <- 1
while(nn < n + 1) {
xcross <- as.matrix(x[-nn,])	
ycross <- as.matrix(y[-nn,])
mca <- svd(t(ycross) %*% xcross)
xtrunc <- xcross %*% mca$v[,1:nmodes]
ytrunc <- ycross %*% mca$u[,1:nmodes]
# Estimate prior parameters yb (1 x q) and c (q x q)
prior <- faprior(ytrunc)
# Estimate likelihood parameters g (p x q), yzerogt (n-1 x p) and
# s (p x q)
lik <- falikelihood(xtrunc, ytrunc)
# Estimate posterior parameters ya (1 x q) and d (q x q)
out <- faprediction(x[nn,] %*% mca$v[,1:nmodes], ytrunc, prior$yb, prior$c, lik$g, lik$yzerogt, lik$s)
ya[nn,] <- out$ya %*% t(mca$u[,1:nmodes])	
d[nn,] <- t(sqrt(diag(mca$u[,1:nmodes] %*% out$d %*% t(mca$u[,1:nmodes]))))
nn <- nn + 1
} # end of cross-validation loop
}
list(ya=ya,d=d)
}