falikelihood <- function(x, y){

# This function estimates the parameters G, yoG' and S
# of the likelihood distribution x|y ~ N([y-yo]G',S), 
# where the symbol ' denotes transpose matrix.
#
# NOTE: in the code all parameters appear
# in lower case letters. 
#
# USAGE: falikelihood(x,y)
#
# 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 var.
#
# OUTPUTS:
# 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 p forecast error covariance matrix
#
# 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 dimentions
#
p <- ncol(x)
q <- ncol(y)
n <- nrow(y)
#
# Define g, yzerogt and s
#
g <- matrix(nrow = p, ncol = q)
yzerogt <- matrix(nrow = n, ncol = p)
s <- matrix(nrow = p, ncol = p)
#
# Calculate column means of x and y
#
# xmean is a 1 x p matrix
#
xmean <- t(apply(x, 2, mean))
#
# ymean is a 1 x q matrix
#
ymean <- t(apply(y, 2, mean))
#
# Estimate covariance matrices
#
Sxx <- var(x)*((n-1)/n)
Syy <- var(y)*((n-1)/n)
Sxy <- var(x, y)*((n-1)/n)
Syx <- t(Sxy)
#
# Inverte Syy
#	
Syyi <- solve(Syy)
#
# Estimate forecast error covariance (s)
#
s <- Sxx - Sxy %*% Syyi %*% Syx
#
# Estimate intercept (yzerogt) and slope (g) from 
# the regression between x and y
#
# g is a p x q matrix
#
g <- Sxy %*% Syyi
#
# yzerogt is a 1 x q matrix. 
#
yzerogt <-  - (xmean - ymean %*% t(g))
list(g = g, yzerogt = yzerogt, s = s)
}