*** Factor Analysis  ***  May 17, 2006.

Monthly log returns of IBM, Hewlett-Packard, Intel, Merrill Lynch, and 
Morgan Stanley Dean Witter from 1990 to 2003.

> setwd("C:/teaching/ama")
> da=read.table("m-pca5c-9003.txt",header=T)
> dim(da)
[1] 168   6
> rtn=da[,1:5]
> mean(rtn)
      IBM       HPQ      INTC       MER       MWD 
0.9710476 1.0353095 2.0366369 1.8732917 1.6031131 

> var(rtn)
          IBM       HPQ      INTC       MER       MWD
IBM  87.28048  48.76290  56.00906  31.00182  38.87554
HPQ  48.76290 137.58537  86.07789  57.93111  60.32278
INTC 56.00906  86.07789 173.61491  52.39257  56.49784
MER  31.00182  57.93111  52.39257 110.30512  93.93719
MWD  38.87554  60.32278  56.49784  93.93719 121.90800

> m1=princomp(rtn)
> summary(m1)
Importance of components:
                           Comp.1     Comp.2    Comp.3     Comp.4     Comp.5
Standard deviation     19.1512557 10.8345557 8.0099130 7.56366645 4.62750533
Proportion of Variance  0.5850173  0.1872390 0.1023363 0.09125127 0.03415608
Cumulative Proportion   0.5850173  0.7722563 0.8745926 0.96584392 1.00000000

> names(m1)
[1] "sdev"  "loadings" "center"  "scale"  "n.obs"  "scores"  "call"    
> m1$loadings

Loadings:
     Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
IBM  -0.300 -0.165 -0.221  0.907 -0.103
HPQ  -0.487 -0.183 -0.765 -0.376       
INTC -0.543 -0.602  0.565 -0.151       
MER  -0.418  0.526  0.146        -0.720
MWD  -0.450  0.547  0.158         0.685

               Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
SS loadings       1.0    1.0    1.0    1.0    1.0
Proportion Var    0.2    0.2    0.2    0.2    0.2
Cumulative Var    0.2    0.4    0.6    0.8    1.0
> L = cbind(19.15*m1$loading[,1],10.83*m1$loading[,2]) % Estimates of 
> L                                                    % factor loading
           [,1]      [,2]
IBM   -5.753797 -1.790197
HPQ   -9.332774 -1.984504
INTC -10.401020 -6.524014
MER   -8.013241  5.698248
MWD   -8.609437  5.925627

> LLT = L%*%t(L)
> LLT
          IBM       HPQ      INTC       MER       MWD
IBM  36.31099  57.25154  71.52463  35.90558  38.92891
HPQ  57.25154  91.03892 110.01729  63.47757  68.59050
INTC 71.52463 110.01729 150.74396  46.17043  50.88805
MER  35.90558  63.47757  46.17043  96.68206 102.75519
MWD  38.92891  68.59050  50.88805 102.75519 109.23546

> psi=diag(var(rtn))-diag(LLT)
> psi
     IBM      HPQ     INTC      MER      MWD 
50.96950 46.54645 22.87094 13.62306 12.67253 

*** Use correlation matrix ****
*** Standardize the return series
> v=var(rtn)
> v1=diag(v)
> v1
      IBM       HPQ      INTC       MER       MWD 
 87.28048 137.58537 173.61491 110.30512 121.90800 
> v1=sqrt(v1)
> v2=diag(v1)
> v2
         [,1]     [,2]    [,3]     [,4]     [,5]
[1,] 9.342402  0.00000  0.0000  0.00000  0.00000
[2,] 0.000000 11.72968  0.0000  0.00000  0.00000
[3,] 0.000000  0.00000 13.1763  0.00000  0.00000
[4,] 0.000000  0.00000  0.0000 10.50262  0.00000
[5,] 0.000000  0.00000  0.0000  0.00000 11.04120

> v2inv=solve(v2)
> srtn=rtn1%*%v2inv
> var(srtn)
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 1.0000000 0.4449844 0.4549946 0.3159590 0.3768789
[2,] 0.4449844 1.0000000 0.5569446 0.4702490 0.4657781
[3,] 0.4549946 0.5569446 1.0000000 0.3785980 0.3883491
[4,] 0.3159590 0.4702490 0.3785980 1.0000000 0.8100720
[5,] 0.3768789 0.4657781 0.3883491 0.8100720 1.0000000
> m2=princomp(srtn)
> summary(m2)
Importance of components:
                          Comp.1    Comp.2    Comp.3     Comp.4     Comp.5
Standard deviation     1.6913269 0.9573868 0.7598213 0.65613597 0.43037104
Proportion of Variance 0.5755432 0.1844156 0.1161571 0.08661847 0.03726567
Cumulative Proportion  0.5755432 0.7599588 0.8761159 0.96273433 1.00000000
> m2$loadings
Loadings:
     Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
[1,]  0.388 -0.471  0.787              
[2,]  0.457 -0.253 -0.413 -0.745       
[3,]  0.425 -0.443 -0.439  0.656       
[4,]  0.474  0.530               -0.699
[5,]  0.484  0.487  0.131         0.709

               Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
SS loadings       1.0    1.0    1.0    1.0    1.0
Proportion Var    0.2    0.2    0.2    0.2    0.2
Cumulative Var    0.2    0.4    0.6    0.8    1.0
> L=cbind(1.691*m2$loadings[,1],.9574*m2$loadings[,2])
> L
          [,1]       [,2]
[1,] 0.6567269 -0.4507447
[2,] 0.7729132 -0.2423579
[3,] 0.7193883 -0.4242094
[4,] 0.8021870  0.5076046
[5,] 0.8183962  0.4659298
> LLT=L%*%t(L)
> LLT
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.6344610 0.6168344 0.6636518 0.2980177 0.3274475
[2,] 0.6168344 0.6561321 0.6588352 0.4969990 0.5196275
[3,] 0.6636518 0.6588352 0.6974731 0.3617533 0.3910929
[4,] 0.2980177 0.4969990 0.3617533 0.9011665 0.8930150
[5,] 0.3274475 0.5196275 0.3910929 0.8930150 0.8868630
> psi=rep(1,5)-diag(LLT)
> psi
[1] 0.36553895 0.34386787 0.30252688 0.09883351 0.11313702

**** MLE
> m2=factanal(rtn,2)   % varimax rotation is used.
> m2

Call:
factanal(x = rtn, factors = 2)

Uniquenesses:
  IBM   HPQ  INTC   MER   MWD 
0.620 0.442 0.449 0.005 0.317 

Loadings:
     Factor1 Factor2
IBM  0.173   0.592  
HPQ  0.310   0.679  
INTC 0.206   0.714  
MER  0.965   0.252  
MWD  0.747   0.354  

               Factor1 Factor2
SS loadings      1.658   1.510
Proportion Var   0.332   0.302
Cumulative Var   0.332   0.633

Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 1.74 on 1 degree of freedom.
The p-value is 0.187 

> m1r=factanal(rtn,2,rotation="none",scores="regression") %store scores
> names(m1r)
 [1] "converged"  "loadings"   "uniquenesses" "correlation"  "criteria"    
 [6] "factors"    "dof"        "method"       "scores"       "STATISTIC"   
[11] "PVAL"       "n.obs"      "call"        

> m1r
Call:
factanal(x = rtn, factors = 2, scores = "regression", rotation = "none")

Uniquenesses:
  IBM   HPQ  INTC   MER   MWD 
0.620 0.442 0.449 0.005 0.317 

Loadings:
     Factor1 Factor2
IBM   0.322   0.526 
HPQ   0.477   0.575 
INTC  0.385   0.635 
MER   0.997         
MWD   0.813   0.146 

               Factor1 Factor2
SS loadings      2.136   1.032
Proportion Var   0.427   0.206
Cumulative Var   0.427   0.633

Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 1.74 on 1 degree of freedom.
The p-value is 0.187 

> dim(m1r$scores)    % the estimated common factors.
[1] 168   2 
> plot(m1r$scores[,1],type='l')
> plot(m1r$scores[,2],type='l')


****** S-Plus analysis
> rtn=pca5c[,1:5]
> m1=factanal(rtn,2)   % principal component method
> summary(m1)
Importance of factors:
                 Factor1   Factor2 
   SS loadings 1.6247936 1.4859191
Proportion Var 0.3249587 0.2971838
Cumulative Var 0.3249587 0.6221425

The degrees of freedom for the model is 1.

Uniquenesses:
       IBM       HPQ      INTC       MER       MWD 
 0.6328239 0.4464351 0.4318128 0.1848429 0.1933727

Loadings:
     Factor1 Factor2 
 IBM 0.208   0.569  
 HPQ 0.315   0.674  
INTC 0.197   0.728  
 MER 0.858   0.280  
 MWD 0.841   0.316  

> m5=factanal(rtn,2,method="mle")  % MLE method
> summary(m5)
Importance of factors:
                 Factor1   Factor2 
   SS loadings 1.6715061 1.4981246
Proportion Var 0.3343012 0.2996249
Cumulative Var 0.3343012 0.6339261

The degrees of freedom for the model is 1.

Uniquenesses:
       IBM      HPQ      INTC          MER       MWD 
 0.6192882 0.442331 0.4491424 0.0002646188 0.3194369

Loadings:
     Factor1 Factor2 
 IBM 0.176   0.591  
 HPQ 0.314   0.678  
INTC 0.210   0.712  
 MER 0.969   0.245  
 MWD 0.747   0.350  

> m6=rotate(m5,rotation="quartimax")   %Some rotation!
> summary(m6)
Importance of factors:
                 Factor1   Factor2 
   SS loadings 1.2997566 1.8698742
Proportion Var 0.2599513 0.3739748
Cumulative Var 0.2599513 0.6339261

The degrees of freedom for the model is 1.

Uniquenesses:
       IBM      HPQ      INTC          MER       MWD 
 0.6192882 0.442331 0.4491424 0.0002646188 0.3194369

Loadings:
     Factor1 Factor2 
 IBM         0.614  
 HPQ 0.177   0.725  
INTC         0.739  
 MER 0.904   0.428  
 MWD 0.666   0.487