faprediction<-function(x, y, yb, c, g, yzerogt, s){

# This function estimates the mean ya and the covariance d
# of the posterior distribution y|x ~ N(ya,d), 
# where ya = yb + [x - yb g' + yo g'] L'
#       d  = (I - Lg) c
#       L' = (gcg' +s)^(-1) gc
#
# The symbol ' indicates transpose matrix
#
# USAGE: faprediction(x, y, yb, c, g, yzerogt, s)
#
# INPUTS:
# 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 vars.
# yb -> 1 x q matrix containing column means of y
# c  -> q x q observable state covariance matrix
# g  -> p x q matrix containing the calibration operator
# yzerogt -> 1 x q matrix containing the intercept of the
#            regression between x and y
# s  -> p x q forecast error covariance matrix
#
# OUTPUTS:
#
# ya -> 1 x q matrix containing q predicted observable state time series
# d  -> q x q matrix containing predicted observable error covaraiance
#
# Authors: Caio Coelho <c.a.d.s.coelho@reading.ac.uk>
#          David Stephenson <d.b.stephenson@reading.ac.uk>
#
       
       
# Find out p, q and n
#
p <- ncol(x)
q <- ncol(y)
n <- nrow(y)
#
# Define matrices ya and d
#
ya <- matrix(nrow = n, ncol = q)
d <- matrix(nrow = q, ncol = q)
#
# Transpose g
#
gt <- t(g)
#
# Calculate gain/weight matrix transpose (Lt)
#
Lt <- solve(g %*% c %*% gt + s) %*% g %*% c
#
# Calculate posterior mean (ya)
#
ya <- yb + (x - yb %*% gt + yzerogt) %*% Lt
#
# Calculate the posterior covariance matrix (d)
#
d <- (diag(q) - t(Lt) %*% g) %*% c
#
list(ya = ya, d = d)
}