*
* Example 9.6.2 from pp 437-440
*
open data m-apca0103.txt
calendar(panelobs=36,m) 2001
data(format=free,org=columns) 1//2001:01 40//2003:12 id date retn
*
compute nper=36
compute nvar=40
*
* Repackage data into an nper x nvar matrix
*
dec rect rmat(nper,nvar)
ewise rmat(i,j)=retn((j-1)*nper+i)
*
* Do Bai-Ng analysis of number of factors
*
@BaiNg(max=10,center) rmat
*
* Pull the means out of each column and compute the eigenvectors of the outer
* product.
*
dec vect sigsq(nvar)
dec rect workmat(nper,nvar)
compute colavg=%fill(1,nper,1.0/nper)*rmat
ewise workmat(i,j)=rmat(i,j)-colavg(1,j)
*
* Do the OLS estimates with six factors
*
compute nf=6
*
compute omega=%outerxx(workmat)
eigen omega eigval eigvect
*
dec rect x(nper,nf+1)
ewise x(i,j)=%if(j==1,1.0,eigvect(i,j-1))
compute mxx=x*inv(tr(x)*x)*tr(x)
*
do j=1,nvar
   compute y=%xcol(workmat,j)
   compute sigsq(j)=(%normsqr(y)-%qform(mxx,y))/nper
   compute lambda=inv(tr(x)*x)*tr(x)*y
end do j
*
* GLS
*
disp "GLS"
*
ewise workmat(i,j)=workmat(i,j)/sqrt(sigsq(j))
compute omega=%outerxx(workmat)
eigen omega eigval eigvect
ewise x(i,j)=%if(j==1,1.0,eigvect(i,j-1))
compute mxx=x*inv(tr(x)*x)*tr(x)
*
do j=1,nvar
   compute y=%xcol(workmat,j)*sqrt(sigsq(j))
   compute lambda=(inv(tr(x)*x)*tr(x)*y)
   compute sigsq(j)=(%normsqr(y)-%qform(mxx,y))/nper
   disp lambda(2,1) 1-nper*sigsq(j)/%normsqr(y)
end do j