# creates a mean forecast value for each year at all stations by dropping each
# year and using the rest of the data to forecast that year
# x is the independent variable, y is the dependent variable (forecast y from x)
## performs PCA of x and then does CCA on the first 15 PCs of x with y
###

CCAForecast=function(x, y)
{

N = dim(x)[1]
M=dim(x)[2]
J=dim(y)[2]
J=min(M,J)
ypred=matrix(0,nrow=N,ncol=J)
index=1:N
for(i in 1:N){
#drop a point..
index1=index[index != i]
xp=x[i,] #prediction point
#get the rest of the data..
X=x[index1,]
xMean=colMeans(X)
xSd=apply(X,2,sd)
Y=y[index1,]
yMean=colMeans(Y)
ySd=apply(Y,2,sd)
## scale the data to be predicted..
xp = (xp-xMean)/xSd
# get the first 6 PCS of X (SST)
X=scale(X)

zs=var(X)
zsvd=svd(zs)
pcs=t(t(zsvd$u) %*% t(X))

#xPca = myPCA(X)
#xpcs=xPca$pcs[,1:6]

xpcs = pcs[,1:15]
X=xpcs
### Perform CCA on the 15 PCS of y (SST) and six streamflows..
X1 = scale(X)
Y1=scale(Y)
Qx1 = qr.Q(qr(X1))
Qy1 = qr.Q(qr(Y1))
VV=t(Qx1) %*% Qy1
T11 = qr.R(qr(X1))
T22 = qr.R(qr(Y1))
BB = solve(T22) %*% svd(VV)$v * sqrt(length(Y1[,1]-1))
BB = svd(VV)$v
wm1=Y1 %*% BB
Fyy = var(Y1) %*% BB
AA = solve(T11) %*% svd(VV)$u * sqrt(length(Y1[,1]-1))
AA = svd(VV)$u
m1 = X1 %*% AA
cancorln = svd(VV)$d[1:J]
#### get the regression coefficient..
betahat = solve(t(AA) %*% t(X1)%*% X1 %*% AA) %*% t(AA) %*% t(X1) %*% Y1
## get the first six Pc values of xp
vmp = (xp %*% zsvd$u)[,1:15]
vmp = vmp %*% AA
### predict
yp = vmp %*% betahat
### scale back
yp = yp*ySd + yMean
ypred[i,]=yp
}

ypred

}