In [1]:
# rm(list=ls())
options(OutDec = ",")
#==============================================================================
# Dados da bolsa de valores de 'New York' em d07_acoes.txt
#==============================================================================
dados <- read.table("../dados/d07_acoes.txt",header=F,sep=" ", # Tab
col.names=c("JP","Citi","Fargo","Royal","Exxon"))
print(head(dados))
JP Citi Fargo Royal Exxon 1 0,0130338 -0,0078431 -0,0031889 -0,0447693 0,0052151 2 0,0084862 0,0166886 -0,0062100 0,0119560 0,0134890 3 -0,0179153 -0,0086393 0,0100360 0,0000000 -0,0061428 4 0,0215589 -0,0034858 0,0174353 -0,0285917 -0,0069534 5 0,0108225 0,0037167 -0,0101345 0,0291900 0,0409751 6 0,0101713 -0,0121978 -0,0083768 0,0137083 0,0029895
In [2]:
#==============================================================================
# Correlação amostral
#==============================================================================
n <- dim(dados)[1]
p <- dim(dados)[2]
R <- cor(dados)
In [3]:
#==============================================================================
# Analise fatorial via maxima verossimilhanca
#==============================================================================
AFMV.1 <- factanal(dados,factors=1,rotation="none")
AFMV.2 <- factanal(dados,factors=2,rotation="none")
In [4]:
#==============================================================================
# Estimacao da matriz de cargas dos fatores
#==============================================================================
L.1 <- AFMV.1$loadings[,1]
L.2 <- AFMV.2$loadings[,1:2]
print(L.1)
print(L.2)
JP Citi Fargo Royal Exxon
0,7154389 0,8750793 0,6630512 0,3460438 0,2792017
Factor1 Factor2
JP 0,1205972 0,754267060
Citi 0,3284924 0,785749642
Fargo 0,1876017 0,650217058
Royal 0,9974724 -0,007103505
Exxon 0,6851746 0,026317443
In [5]:
#==============================================================================
# Estimacao das comunalidades
#==============================================================================
sum2 <- function(x){return(sum(x^2))}
h.1 <- L.1^2 # um fator
h.2 <- apply(L.2,1,"sum2") # dois fatores
print(h.1)
print(h.2)
JP Citi Fargo Royal Exxon
0,51185287 0,76576377 0,43963688 0,11974629 0,07795361
JP Citi Fargo Royal Exxon
0,5834625 0,7253097 0,4579766 0,9950016 0,4701569
In [6]:
#==============================================================================
# Estimacao das variancias especificas
#==============================================================================
psi.1 <- AFMV.1$uniqueness
psi.2 <- AFMV.2$uniqueness
print(psi.1)
print(psi.2)
JP Citi Fargo Royal Exxon
0,4881494 0,2342357 0,5603556 0,8802602 0,9220490
JP Citi Fargo Royal Exxon
0,4165374 0,2746902 0,5420233 0,0050000 0,5298429
In [7]:
#==============================================================================
# Estimacao/Recomposicao de R para m=2
#==============================================================================
R.est.1 <- L.1%*%t(L.1)+diag(psi.1)
R.est.2 <- L.2%*%t(L.2)+diag(psi.2)
print(R.est.1)
print(R.est.2)
JP Citi Fargo Royal Exxon
[1,] 1,0000022 0,6260658 0,4743726 0,24757319 0,19975179
[2,] 0,6260658 0,9999994 0,5802224 0,30281574 0,24432365
[3,] 0,4743726 0,5802224 0,9999925 0,22944473 0,18512504
[4,] 0,2475732 0,3028157 0,2294447 1,00000651 0,09661602
[5,] 0,1997518 0,2443237 0,1851250 0,09661602 1,00000263
JP Citi Fargo Royal Exxon
JP 0,9999999 0,6322803 0,5130616 0,1149345 0,1024805
Citi 0,6322803 1,0000000 0,5725336 0,3220805 0,2457536
Fargo 0,5130616 0,5725336 0,9999999 0,1825087 0,1456520
Royal 0,1149345 0,3220805 0,1825087 1,0000016 0,6832558
Exxon 0,1024805 0,2457536 0,1456520 0,6832558 0,9999997
In [8]:
#==============================================================================
# Matriz residual
#==============================================================================
U <- R - R.est.2
print(U)
JP Citi Fargo Royal Exxon JP 1,055860e-07 7,496780e-06 -2,564223e-03 -3,325561e-04 5,198222e-02 Citi 7,496780e-06 3,255672e-08 1,608871e-03 2,116218e-04 -3,307885e-02 Fargo -2,564223e-03 1,608871e-03 5,157367e-08 -9,518792e-06 5,547153e-04 Royal -3,325561e-04 2,116218e-04 -9,518792e-06 -1,559500e-06 1,218853e-04 Exxon 5,198222e-02 -3,307885e-02 5,547153e-04 1,218853e-04 2,670491e-07
In [9]:
#==============================================================================
# Comparacoes: note as entradas [1,3] das matrizes
#==============================================================================
print(R)
print(R.est.1)
print(R.est.2)
JP Citi Fargo Royal Exxon
JP 1,0000000 0,6322878 0,5104973 0,1146019 0,1544628
Citi 0,6322878 1,0000000 0,5741424 0,3222921 0,2126747
Fargo 0,5104973 0,5741424 1,0000000 0,1824992 0,1462067
Royal 0,1146019 0,3222921 0,1824992 1,0000000 0,6833777
Exxon 0,1544628 0,2126747 0,1462067 0,6833777 1,0000000
JP Citi Fargo Royal Exxon
[1,] 1,0000022 0,6260658 0,4743726 0,24757319 0,19975179
[2,] 0,6260658 0,9999994 0,5802224 0,30281574 0,24432365
[3,] 0,4743726 0,5802224 0,9999925 0,22944473 0,18512504
[4,] 0,2475732 0,3028157 0,2294447 1,00000651 0,09661602
[5,] 0,1997518 0,2443237 0,1851250 0,09661602 1,00000263
JP Citi Fargo Royal Exxon
JP 0,9999999 0,6322803 0,5130616 0,1149345 0,1024805
Citi 0,6322803 1,0000000 0,5725336 0,3220805 0,2457536
Fargo 0,5130616 0,5725336 0,9999999 0,1825087 0,1456520
Royal 0,1149345 0,3220805 0,1825087 1,0000016 0,6832558
Exxon 0,1024805 0,2457536 0,1456520 0,6832558 0,9999997
In [10]:
#==============================================================================
# Soma de quadrado das cargas para o k-esimo fator
#==============================================================================
ssc.1 <- sum2(L.1) # um fator
ssc.2 <- apply(L.2,2,"sum2") # dois fatores
print(ssc.1)
print(ssc.2)
[1] 1,914953 Factor1 Factor2 1,622061 1,609847
In [11]:
#==============================================================================
# Proporcao da variabilidade total
#==============================================================================
print(ssc.1/p) # um fator
print(ssc.2/p) # dois fatores
prop.ac <- cumsum(ssc.2)/p
print("==== # Proporção acumulada para dois fatores ====")
print(prop.ac)
[1] 0,3829907 Factor1 Factor2 0,3244121 0,3219693 [1] "==== # Proporção acumulada para dois fatores ====" Factor1 Factor2 0,3244121 0,6463815
In [12]:
#==============================================================================
# Para comparar os resultados acima, podemos ver diretamente:
#==============================================================================
print(AFMV.1)
print(AFMV.2)
Call:
factanal(x = dados, factors = 1, rotation = "none")
Uniquenesses:
JP Citi Fargo Royal Exxon
0,488 0,234 0,560 0,880 0,922
Loadings:
Factor1
JP 0,715
Citi 0,875
Fargo 0,663
Royal 0,346
Exxon 0,279
Factor1
SS loadings 1,915
Proportion Var 0,383
Test of the hypothesis that 1 factor is sufficient.
The chi square statistic is 62,22 on 5 degrees of freedom.
The p-value is 4,22e-12
Call:
factanal(x = dados, factors = 2, rotation = "none")
Uniquenesses:
JP Citi Fargo Royal Exxon
0,417 0,275 0,542 0,005 0,530
Loadings:
Factor1 Factor2
JP 0,121 0,754
Citi 0,328 0,786
Fargo 0,188 0,650
Royal 0,997
Exxon 0,685
Factor1 Factor2
SS loadings 1,622 1,610
Proportion Var 0,324 0,322
Cumulative Var 0,324 0,646
Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 1,97 on 1 degree of freedom.
The p-value is 0,16
In [13]:
#==============================================================================
# Testes de adequacao
#==============================================================================
nu.1 <- ((p-1)^2-p-1)/2
nu.2 <- ((p-2)^2-p-2)/2
print(nu.1)
print(nu.2)
[1] 5 [1] 1
In [14]:
#==============================================================================
# Estatística do teste (1)
#==============================================================================
m1 <- 1
W1 <- (n-1-(2*p + 4*m1 + 5)/6)*( log(det(R.est.1))-log(det(R)))
print(W1)
[1] 62,22057
In [15]:
#==============================================================================
# Resultado do teste (1)
#==============================================================================
print(W1 > qchisq(0.95,nu.1))
valorp.1 <- 1-pchisq(W1,nu.1)
print(valorp.1)
[1] TRUE [1] 4,221401e-12
In [16]:
#==============================================================================
# Compare com o resultado de
#==============================================================================
print(AFMV.1)
Call:
factanal(x = dados, factors = 1, rotation = "none")
Uniquenesses:
JP Citi Fargo Royal Exxon
0,488 0,234 0,560 0,880 0,922
Loadings:
Factor1
JP 0,715
Citi 0,875
Fargo 0,663
Royal 0,346
Exxon 0,279
Factor1
SS loadings 1,915
Proportion Var 0,383
Test of the hypothesis that 1 factor is sufficient.
The chi square statistic is 62,22 on 5 degrees of freedom.
The p-value is 4,22e-12
In [17]:
#==============================================================================
# Estatística do teste (2)
#==============================================================================
m2 <- 2
W2 <- (n-1-(2*p + 4*m2 + 5)/6)*( log(det(R.est.2))-log(det(R)))
print(W2)
[1] 2,004672
In [18]:
#==============================================================================
# Resultado do teste (2)
#==============================================================================
print(W2 > qchisq(0.95,nu.2))
valorp.2 <- 1-pchisq(W2,nu.2)
print(valorp.2)
[1] FALSE [1] 0,1568152
In [19]:
#==============================================================================
# Compare com o resultado de
#==============================================================================
print(AFMV.2)
Call:
factanal(x = dados, factors = 2, rotation = "none")
Uniquenesses:
JP Citi Fargo Royal Exxon
0,417 0,275 0,542 0,005 0,530
Loadings:
Factor1 Factor2
JP 0,121 0,754
Citi 0,328 0,786
Fargo 0,188 0,650
Royal 0,997
Exxon 0,685
Factor1 Factor2
SS loadings 1,622 1,610
Proportion Var 0,324 0,322
Cumulative Var 0,324 0,646
Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 1,97 on 1 degree of freedom.
The p-value is 0,16
In [20]:
#==============================================================================
# Rotacao de fatores
#==============================================================================
# Ver esta parte so' apos estudar rotacao de fatores
# Poderiamos utilizar a funcao "varimax" em L.2
# Mas utilizaremos a opcao "varimax" em "factanal"
#==============================================================================
AFMVmax.2 <- factanal(dados,factors=2,rotation="varimax")
In [21]:
#==============================================================================
# Estimacao da matriz de cargas dos fatores
#==============================================================================
Lmax.2 <- AFMVmax.2$loadings[,1:2]
print(Lmax.2)
Factor1 Factor2 JP 0,7632891 0,02919252 Citi 0,8194973 0,23180593 Fargo 0,6680337 0,10820147 Royal 0,1126721 0,99111380 Exxon 0,1083670 0,67706236
In [22]:
#==============================================================================
# Estimacao dos escores dos fatores
#==============================================================================
AFMV.2.Reg <- factanal(dados,factors=2,rotation="none",scores="regression")
AFMV.2.Bar <- factanal(dados,factors=2,rotation="none",scores="Bartlett")
par(mfrow=c(1,1),lab=c(9,9,5))
plot(AFMV.2.Reg$scores[,1],AFMV.2.Bar$scores[,1],pch=15,bty="n",
xlim=c(-4,4),ylim=c(-4,4),
xlab=expression(paste(hat(f)[1],(Reg),sep="")),
ylab=expression(paste(hat(f)[1],(Bar),sep="")))
par(mfrow=c(1,1),lab=c(9,9,5))
plot(AFMV.2.Reg$scores[,2],AFMV.2.Bar$scores[,2],pch=15,bty="n",
xlim=c(-4,4),ylim=c(-4,4),
xlab=expression(paste(hat(f)[2],(Reg),sep="")),
ylab=expression(paste(hat(f)[2],(Bar),sep="")))
In [23]:
#==============================================================================
# Correlacao
#==============================================================================
print(cor(AFMV.2.Reg$scores,AFMV.2.Bar$scores))
Factor1 Factor2 Factor1 1,000000e+00 1,699021e-17 Factor2 -8,416802e-18 1,000000e+00
In [24]:
#==============================================================================
# Verificando a estabilidade da solucao
# Dados completos e subconjuntos de 51 e 52 observacoes
#==============================================================================
AFMVmax.2a <- factanal(dados[1:52,],factors=2,rotation="varimax")
AFMVmax.2b <- factanal(dados[53:n,],factors=2,rotation="varimax")
Lmax.2a <- AFMVmax.2a$loadings[,1:2]
Lmax.2b <- AFMVmax.2b$loadings[,1:2]
print("==== Solução com todos os dados ====")
print(AFMVmax.2)
print("==== Solução com a primeira parte dos dados ====")
print(AFMVmax.2a)
print("==== Solução com a segunda parte dos dados ====")
print(AFMVmax.2b)
[1] "==== Solução com todos os dados ===="
Call:
factanal(x = dados, factors = 2, rotation = "varimax")
Uniquenesses:
JP Citi Fargo Royal Exxon
0,417 0,275 0,542 0,005 0,530
Loadings:
Factor1 Factor2
JP 0,763
Citi 0,819 0,232
Fargo 0,668 0,108
Royal 0,113 0,991
Exxon 0,108 0,677
Factor1 Factor2
SS loadings 1,725 1,507
Proportion Var 0,345 0,301
Cumulative Var 0,345 0,646
Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 1,97 on 1 degree of freedom.
The p-value is 0,16
[1] "==== Solução com a primeira parte dos dados ===="
Call:
factanal(x = dados[1:52, ], factors = 2, rotation = "varimax")
Uniquenesses:
JP Citi Fargo Royal Exxon
0,362 0,273 0,747 0,005 0,496
Loadings:
Factor1 Factor2
JP 0,799
Citi 0,813 0,258
Fargo 0,469 0,183
Royal 0,119 0,990
Exxon 0,180 0,687
Factor1 Factor2
SS loadings 1,565 1,553
Proportion Var 0,313 0,311
Cumulative Var 0,313 0,624
Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 0,62 on 1 degree of freedom.
The p-value is 0,43
[1] "==== Solução com a segunda parte dos dados ===="
Call:
factanal(x = dados[53:n, ], factors = 2, rotation = "varimax")
Uniquenesses:
JP Citi Fargo Royal Exxon
0,447 0,299 0,156 0,005 0,526
Loadings:
Factor1 Factor2
JP 0,741
Citi 0,805 0,231
Fargo 0,918
Royal 0,110 0,991
Exxon 0,684
Factor1 Factor2
SS loadings 2,058 1,510
Proportion Var 0,412 0,302
Cumulative Var 0,412 0,714
Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 1,33 on 1 degree of freedom.
The p-value is 0,248