/* first input matrix XY */
scalar n = rowsof(XY)
mat u = J(n,1,1)
mat mv = (1/n)*u'*XY
mat mm = u*mv
mat d = XY-mm

/* covariance matrix */
mat s=(1/(n-1))*d'*d    

/* partition covariance matrix */
mat syy = s[1..3,1..3]
mat sxx = s[4..8,4..8]
mat syx = s[1..3,4..8]
mat sxy = syx'

mat t1 = syminv(syy)*syx
mat t2 = syminv(sxx)*sxy

mat C = syx*syminv(sxx)*sxy
mat F = cholesky(syminv(syy))
mat D = F'*C*F
mat symeigen W L = D
mat U = F*W

mat cc = cholesky(diag(L)) /* canonical correlations */
mat ccc = syminv(cc)
mat cc = vecdiag(cc)       /* make row vector */
mat V = t2*U*ccc


mat list L, title(eigenvalues)
mat list cc, title(canonical correlations)
mat list U, title(eigenvectors1)
mat list V, title(eigenvectors2)