
function res= canonical(F,G);
% F is the latent factors 
% G is the observed factors
% p is the order of the deterministic terms (0=constant, 1=constant+trend)


[T,k]=size(F);
[T,m]=size(G);
reg=ones(T,1);
se=zeros(m,1);
se_k=zeros(m,1);
r2=zeros(m,1);
kappa=zeros(m,1);
MLE=zeros(m,1);


 s00=F'*F/T;
 s11=G'*G/T;
 s10=G'*F/T;
sigma=inv(s11)*s10*inv(s00)*s10';


% canonical correlations are the square-roots of the first r eigenvalues 
% of inv(s11)s10 inv(s00)s10'
%[a,b,c]=svd(sigma);
s=eig(sigma);  
s=sort(s,1);  % smallest to largest
r2=flipud(s);   % these are eigenvalues from largest to smallest


f=demean(F);
% compute confidence interval
for i=1:min([k;m]);
  kappa(i)= mean(f(:,i).^4)/(mean(f(:,i).^2)^2)-3;
%  se(i)=1.96*2*sqrt(r2(i)*(1-r2(i)^2))/sqrt(T);
  se(i)=1.96*2*sqrt(r2(i))*(1-r2(i))/sqrt(T);
  se_k(i)=(1+kappa(i)/3)*se(i);

  end;


%res.r2=r2(min([k;m]));
%res.se=se(min([k;m]));
%res.sek=se_k(min([k;m]));

res.r2=r2;
res.se=se;
res.se_k=se_k;



%MLE=-(T/2)*repmat(log(det(s00)),m,1);
%for i=1:m;
%  MLE(i)=MLE(i)-T/2*sum(log(1-r2(i:m)));
%  end;
  
%[a,test_k]=max(MLE);





 
