test=matrix(scan("Leesferry-mon-data.txt"),ncol=13,byrow=T)
test=test[,2:13]

xt1=test[,1:2]		#Jan-Feb

xt=test[,3:4]		#Mar-Apr
nyrs=length(xt[,1])
nyrs1=nyrs-1

xt=scale(xt,scale=FALSE)
xt1=scale(xt1,scale=FALSE)


Mo = var(xt)      #Sxx
Mo = var(xt)      #Sxx
Moj1 = var(xt1)      #Sxt1

#OR

Mo = (1/nyrs1)*xt%*%t(xt)
Moj1 = (1/nyrs1)*xt1%*%t(xt1)


M1 = var(xt, xt1)   #Syx

#OR

M1 = (1/nyrs1)*t(xt)%*%xt1

A= M1%*%solve(Moj1) # Syx%*%solve(Sxx)

D = Mo - (M1%*%solve(Moj1))%*%(t(M1))
zz=svd(D)
B=zz$u %*% diag(sqrt(zz$d)) 	#gives a m x m matrix
# zz$u is the eigenvector and diag(sqrt(zz$d)) are the eigen values

W=matrix(rnorm(1,0,1),nrow=2,ncol=1)

xsim=c(mean(xt1[,1],xt1[,2])
xsim=A %*% t(xsim) + B %*% W