# 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
#==============================================================================
# Correlação amostral
#==============================================================================
n <- dim(dados)[1]
p <- dim(dados)[2]
R <- cor(dados)
#==============================================================================
# Analise fatorial via maxima verossimilhanca
#==============================================================================
AFMV.1 <- factanal(dados,factors=1,rotation="none")
AFMV.2 <- factanal(dados,factors=2,rotation="none")
#==============================================================================
# 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
#==============================================================================
# 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
#==============================================================================
# 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
#==============================================================================
# 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
#==============================================================================
# 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,255673e-08 1,608871e-03 2,116218e-04 -3,307885e-02 Fargo -2,564223e-03 1,608871e-03 5,157370e-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
#==============================================================================
# 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
#==============================================================================
# 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
#==============================================================================
# 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
#==============================================================================
# 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
#==============================================================================
# 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
#==============================================================================
# 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
#==============================================================================
# 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
#==============================================================================
# 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
#==============================================================================
# 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
#==============================================================================
# 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
#==============================================================================
# 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