/* 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)