Vamos considerar os exemplos do professor em sala de aula e estudar um pouco sobre o tópico de análise fatorial.
\[\Sigma = \begin{pmatrix} 19 & 30 & 2 & 12 \\ 30 & 57 & 5 & 23 \\ 2 & 5 & 38 & 47\\ 12 & 23 & 47 & 68 \end{pmatrix}\]
verifique que \(\Sigma = LL^T + \Psi\) para m = 2, com
\[L = \begin{pmatrix} 4 & 1 \\ 7 & 2 \\ -1 & 6\\ 1 & 8 \end{pmatrix}\]
e \(\Psi = diag(2,4,1,3).\)
SIGMA = matrix(c(19,30,2,12,30,57,5,23,2,5,38,47,12,23,47,68),4,4,byrow = T)
L = matrix(c(4,1,7,2,-1,6,1,8),4,2,byrow = T)
PSI = diag(c(2,4,1,3))
SIGMA
## [,1] [,2] [,3] [,4]
## [1,] 19 30 2 12
## [2,] 30 57 5 23
## [3,] 2 5 38 47
## [4,] 12 23 47 68
L
## [,1] [,2]
## [1,] 4 1
## [2,] 7 2
## [3,] -1 6
## [4,] 1 8
PSI
## [,1] [,2] [,3] [,4]
## [1,] 2 0 0 0
## [2,] 0 4 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 3
Vamos verificar, se de fato essas cargas fatoriais, com esses dois fatores estão adequados:
L%*%t(L)+PSI # Verificando os valores
## [,1] [,2] [,3] [,4]
## [1,] 19 30 2 12
## [2,] 30 57 5 23
## [3,] 2 5 38 47
## [4,] 12 23 47 68
Vamos comparar se é igual a Matriz de covariâncias Sigma.
SIGMA - (L%*%t(L) + PSI)
## [,1] [,2] [,3] [,4]
## [1,] 0 0 0 0
## [2,] 0 0 0 0
## [3,] 0 0 0 0
## [4,] 0 0 0 0
Vemos, que de fato isto foi atendido.
comunal = function(x) { sum(x^2) }
h = apply(L, 1, "comunal")
h
## [1] 17 53 37 65
h+diag(PSI)
## [1] 19 57 38 68
h+diag(PSI)==diag(SIGMA)
## [1] TRUE TRUE TRUE TRUE
De fato, a soma da comunalidade com a variância específica é igual a variância que se encontra na diagonal da matriz de covariância.
\[\Sigma = \begin{pmatrix} 1 & 0.9 & 0.7 \\ 0.9 & 1 & 0.4 \\ 0.7 & 0.4 & 1 \end{pmatrix}\]
# A matriz RHO é positiva definida
RHO <- matrix(c(1,0.9,0.7,0.9,1,0.4,0.7,0.4,1),3,3,byrow = T)
RHO
## [,1] [,2] [,3]
## [1,] 1.0 0.9 0.7
## [2,] 0.9 1.0 0.4
## [3,] 0.7 0.4 1.0
Pode ser obtida uma solução através de componentes principais, todavia ela não é exata, só é exata se tivermos m=p=3, que não haveria necessidade já que o objetivo é diminuir o número de dimensões.
dec <- eigen(RHO)
delta <- dec$values #autovalores
P <- dec$vectors #autovetores
ell.1 <- sqrt(delta[1])*P[,1]
ell.2 <- sqrt(delta[2])*P[,2]
ell.3 <- sqrt(delta[3])*P[,3]
L <- cbind(ell.1,ell.2,ell.3)
L
## ell.1 ell.2 ell.3
## [1,] -0.9874197 -0.0871829 0.13191463
## [2,] -0.8846479 -0.4553643 -0.10020701
## [3,] -0.7720339 0.6332923 -0.05389295
RHO-ell.1%*%t(ell.1) # Se considerar m=1 a solucao nao é exata
## [,1] [,2] [,3]
## [1,] 0.02500233 0.02648121 -0.06232153
## [2,] 0.02648121 0.21739805 -0.28297823
## [3,] -0.06232153 -0.28297823 0.40396359
# Aproveitamos o exemplo para verificar:
delta/3 # Proporcao da variabilidade
## [1] 0.78454534 0.20533887 0.01011579
cumsum(delta)/3 # Proporcao da variabilidade acumulada
## [1] 0.7845453 0.9898842 1.0000000
# Para conferir, note que delta_k = ell.1_k'*ell.1_k
# Logo, podemos ver a proporcao da variabilidade devida ao
# k-esimo fator atraves da soma de quadrados dos
# valores da k-esima coluna de L
delta
## [1] 2.35363603 0.61601660 0.03034736
cbind((ell.1)%*%ell.1,(ell.2)%*%ell.2,(ell.3)%*%ell.3)
## [,1] [,2] [,3]
## [1,] 2.353636 0.6160166 0.03034736
\[R = \begin{pmatrix} 1 & 0.02 & 0.96 & 0.42 & 0.01 \\ 0.02 & 1 & 0.13 & 0.71 & 0.85 \\ 0.96 & 0.13 & 1 & 0.50 & 0.11 \\ 0.42 & 0.71 & 0.50 & 1 & 0.79 \\ 0.01 & 0.85 & 0.11 & 0.79 & 1 \end{pmatrix}\]
Pode-se observar que as variáveis 1 e 3 (gosto e aroma) formam um grupo, assim como as variáveis 2 e 5 (preço e energia). A variável 4, “bom para lanche” está mais próxima do grupo (preço, energia).
R <- matrix(c(
1, 0.02, 0.96, 0.42, 0.01,
0.02, 1, 0.13, 0.71, 0.85,
0.96, 0.13, 1, 0.50, 0.11,
0.42, 0.71, 0.50, 1, 0.79,
0.01, 0.85, 0.11, 0.79, 1),5,5, byrow = T)
R
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.00 0.02 0.96 0.42 0.01
## [2,] 0.02 1.00 0.13 0.71 0.85
## [3,] 0.96 0.13 1.00 0.50 0.11
## [4,] 0.42 0.71 0.50 1.00 0.79
## [5,] 0.01 0.85 0.11 0.79 1.00
p <- 5 # número de variáveis
AF <- eigen(R)
delta <- AF$values #autovalores
P <- AF$vectors #autovetores
delta
## [1] 2.85309042 1.80633245 0.20449022 0.10240947 0.03367744
round(P,3)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.331 -0.607 0.098 0.139 0.702
## [2,] 0.460 0.390 0.743 -0.282 0.072
## [3,] 0.382 -0.557 0.168 0.117 -0.709
## [4,] 0.556 0.078 -0.602 -0.568 0.002
## [5,] 0.473 0.404 -0.221 0.751 0.009
delta/p # Proporcao da variabilidade
## [1] 0.570618084 0.361266491 0.040898045 0.020481893 0.006735487
Verificando a proporção da variabilidade explicada para saber quantos fatores vamos utilizar.
cumsum(delta)/p
## [1] 0.5706181 0.9318846 0.9727826 0.9932645 1.0000000
Como podemos enxergar, com dois fatores temos uma explicação de cerca de 93.19% da variabilidade dos dados. Sendo assim, vamos escolher m = 2.
ell.1 <- sqrt(delta[1])*P[,1]
ell.2 <- sqrt(delta[2])*P[,2]
L <- cbind(ell.1,ell.2)
L
## ell.1 ell.2
## [1,] 0.5598618 -0.8160981
## [2,] 0.7772594 0.5242021
## [3,] 0.6453364 -0.7479464
## [4,] 0.9391057 0.1049187
## [5,] 0.7982069 0.5432281
L%*%t(L)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.979461354 0.007357539 0.9716968 0.4401455 0.003558179
## [2,] 0.007357539 0.878920022 0.1095187 0.7849273 0.905175175
## [3,] 0.971696840 0.109518709 0.9758829 0.5275656 0.108806485
## [4,] 0.440145533 0.784927343 0.5275656 0.8929275 0.806595488
## [5,] 0.003558179 0.905175175 0.1088065 0.8065955 0.932231117
round(R-L%*%t(L),digits = 4)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.0205 0.0126 -0.0117 -0.0201 0.0064
## [2,] 0.0126 0.1211 0.0205 -0.0749 -0.0552
## [3,] -0.0117 0.0205 0.0241 -0.0276 0.0012
## [4,] -0.0201 -0.0749 -0.0276 0.1071 -0.0166
## [5,] 0.0064 -0.0552 0.0012 -0.0166 0.0678
# Para conferir, note que delta_j = ell.1_j'*ell.1_j
# Logo, podemos ver a proporcao da variabilidade devida ao
# k-esimo fator atraves da soma de quadrados dos
# valores da k-esima coluna de L
delta
## [1] 2.85309042 1.80633245 0.20449022 0.10240947 0.03367744
cbind((ell.1)%*%ell.1,(ell.2)%*%ell.2)
## [,1] [,2]
## [1,] 2.85309 1.806332
# Estimacao das comunalidades
sum2 <- function(x){return(sum(x^2))}
h <- apply(L,1,"sum2")
h
## [1] 0.9794614 0.8789200 0.9758829 0.8929275 0.9322311
# Estimacao das variancias especificas
psi <- 1 - h # variancia das variaveis padronizadas - comundalidades
psi
## [1] 0.02053865 0.12107998 0.02411712 0.10707250 0.06776888
R_estimado <- L%*%t(L) + diag(psi)
R_estimado
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.000000000 0.007357539 0.9716968 0.4401455 0.003558179
## [2,] 0.007357539 1.000000000 0.1095187 0.7849273 0.905175175
## [3,] 0.971696840 0.109518709 1.0000000 0.5275656 0.108806485
## [4,] 0.440145533 0.784927343 0.5275656 1.0000000 0.806595488
## [5,] 0.003558179 0.905175175 0.1088065 0.8065955 1.000000000
#matriz residual
U <- R - R_estimado
U
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.000000000 0.01264246 -0.011696840 -0.02014553 0.006441821
## [2,] 0.012642461 0.00000000 0.020481291 -0.07492734 -0.055175175
## [3,] -0.011696840 0.02048129 0.000000000 -0.02756557 0.001193515
## [4,] -0.020145533 -0.07492734 -0.027565574 0.00000000 -0.016595488
## [5,] 0.006441821 -0.05517518 0.001193515 -0.01659549 0.000000000
Obtivemos boa explicação através de dois fatores
V.L <- varimax(L)
varimax(V.L$loadings)
## $loadings
##
## Loadings:
## ell.1 ell.2
## [1,] -0.989
## [2,] 0.937
## [3,] 0.129 -0.979
## [4,] 0.843 -0.427
## [5,] 0.965
##
## ell.1 ell.2
## SS loadings 2.538 2.121
## Proportion Var 0.508 0.424
## Cumulative Var 0.508 0.932
##
## $rotmat
## [,1] [,2]
## [1,] 0.9999997812 -0.0006614858
## [2,] 0.0006614858 0.9999997812
V.L$loadings
##
## Loadings:
## ell.1 ell.2
## [1,] -0.989
## [2,] 0.937
## [3,] 0.130 -0.979
## [4,] 0.843 -0.427
## [5,] 0.965
##
## ell.1 ell.2
## SS loadings 2.539 2.121
## Proportion Var 0.508 0.424
## Cumulative Var 0.508 0.932
Tivemos essas cargas fatoriais e podemos dizer que o primeiro é uma média ponderada, e o segundo também, identificando dois grupos.
cbind(L,V.L$loadings[,1:2])
## ell.1 ell.2 ell.1 ell.2
## [1,] 0.5598618 -0.8160981 0.02125556 -0.98944912
## [2,] 0.7772594 0.5242021 0.93742129 0.01270189
## [3,] 0.6453364 -0.7479464 0.13009861 -0.97926362
## [4,] 0.9391057 0.1049187 0.84310820 -0.42672715
## [5,] 0.7982069 0.5432281 0.96536898 0.01714214
rm(list=ls(all=T))
setwd("C:/Users/Pessoal/Desktop/ESTATÍSTICA/UFPB/8º PERÍODO/ANÁLISE MULTIVARIADA II/AULAS/ANÁLISE FATORIAL")
dados <- read.table("acp-ex4-data.txt", sep = "\t",
dec = ".", header = F,
col.names = c("JP", "Citi", "Fargo", "Royal", "Exxon"))
DescTools::Desc(dados, plotit = F, digits = 4)
## ------------------------------------------------------------------------------
## Describe dados (data.frame):
##
## data frame: 103 obs. of 5 variables
## 103 complete cases (100.0%)
##
## Nr ColName Class NAs Levels
## 1 JP numeric .
## 2 Citi numeric .
## 3 Fargo numeric .
## 4 Royal numeric .
## 5 Exxon numeric .
##
##
## ------------------------------------------------------------------------------
## 1 - JP (numeric)
##
## length n NAs unique 0s mean meanCI'
## 103 103 0 = n 0 0.0011 -0.0030
## 100.0% 0.0% 0.0% 0.0051
##
## .05 .10 .25 median .75 .90 .95
## -0.0330 -0.0265 -0.0136 0.0034 0.0168 0.0238 0.0302
##
## range sd vcoef mad IQR skew kurt
## 0.0943 0.0208 19.5855 0.0217 0.0304 -0.1128 -0.5367
##
## lowest : -0.0459, -0.0447, -0.0374, -0.0345, -0.0331
## highest: 0.0308, 0.0373, 0.0466, 0.0485, 0.0485
##
## ' 95%-CI (classic)
##
## ------------------------------------------------------------------------------
## 2 - Citi (numeric)
##
## length n NAs unique 0s mean meanCI'
## 103 103 0 = n 0 0.0007 -0.0034
## 100.0% 0.0% 0.0% 0.0047
##
## .05 .10 .25 median .75 .90 .95
## -0.0311 -0.0227 -0.0132 0.0017 0.0140 0.0239 0.0394
##
## range sd vcoef mad IQR skew kurt
## 0.1123 0.0209 31.9575 0.0207 0.0273 0.0893 0.3767
##
## lowest : -0.0598, -0.0497, -0.0431, -0.0408, -0.0313
## highest: 0.0411, 0.0441, 0.0457, 0.0523, 0.0525
##
## ' 95%-CI (classic)
##
## ------------------------------------------------------------------------------
## 3 - Fargo (numeric)
##
## length n NAs unique 0s mean meanCI'
## 103 103 0 102 2 0.0016 -0.0013
## 100.0% 0.0% 1.9% 0.0046
##
## .05 .10 .25 median .75 .90 .95
## -0.0187 -0.0164 -0.0081 0.0003 0.0100 0.0206 0.0289
##
## range sd vcoef mad IQR skew kurt
## 0.0769 0.0150 9.2035 0.0130 0.0181 0.3222 0.0108
##
## lowest : -0.0362, -0.0304, -0.0238, -0.0238, -0.0220
## highest: 0.0318, 0.0344, 0.0350, 0.0382, 0.0407
##
## ' 95%-CI (classic)
##
## ------------------------------------------------------------------------------
## 4 - Royal (numeric)
##
## length n NAs unique 0s mean meanCI'
## 103 103 0 = n 1 0.0040 -0.0012
## 100.0% 0.0% 1.0% 0.0093
##
## .05 .10 .25 median .75 .90 .95
## -0.0446 -0.0305 -0.0145 0.0063 0.0222 0.0363 0.0445
##
## range sd vcoef mad IQR skew kurt
## 0.1159 0.0269 6.6383 0.0245 0.0367 -0.1267 -0.4341
##
## lowest : -0.0539, -0.0519, -0.0508, -0.0498, -0.0482
## highest: 0.0510, 0.0569, 0.0582, 0.0616, 0.0620
##
## ' 95%-CI (classic)
##
## ------------------------------------------------------------------------------
## 5 - Exxon (numeric)
##
## length n NAs unique 0s mean meanCI'
## 103 103 0 = n 0 0.0040 -0.0014
## 100.0% 0.0% 0.0% 0.0094
##
## .05 .10 .25 median .75 .90 .95
## -0.0445 -0.0264 -0.0125 0.0052 0.0216 0.0381 0.0509
##
## range sd vcoef mad IQR skew kurt
## 0.1420 0.0277 6.8515 0.0249 0.0342 -0.1173 0.3017
##
## lowest : -0.0636, -0.0629, -0.0620, -0.0583, -0.0557
## highest: 0.0561, 0.0574, 0.0587, 0.0649, 0.0784
##
## ' 95%-CI (classic)
n <- dim(dados)[1] # quantidade de observações
p <- dim(dados)[2] # quantidade de variáveis
R <- cor(dados) # matriz de correlação das variáveis
round(R,4)
## JP Citi Fargo Royal Exxon
## JP 1.0000 0.6323 0.5105 0.1146 0.1545
## Citi 0.6323 1.0000 0.5741 0.3223 0.2127
## Fargo 0.5105 0.5741 1.0000 0.1825 0.1462
## Royal 0.1146 0.3223 0.1825 1.0000 0.6834
## Exxon 0.1545 0.2127 0.1462 0.6834 1.0000
AFMV.1 <- factanal(dados,factors=1,rotation="none") # Com um único fator
AFMV.2 <- factanal(dados,factors=2,rotation="none") # Com dois fatores
L.1 <- AFMV.1$loadings[,1]
L.2 <- AFMV.2$loadings[,1:2]
L.1
## JP Citi Fargo Royal Exxon
## 0.7154389 0.8750793 0.6630512 0.3460438 0.2792017
L.2
## 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
sum2 <- function(x){return(sum(x^2))}
h.1 <- L.1^2 # um fator
h.2 <- apply(L.2,1,"sum2") # dois fatores
h.1
## JP Citi Fargo Royal Exxon
## 0.51185287 0.76576377 0.43963688 0.11974629 0.07795361
h.2
## JP Citi Fargo Royal Exxon
## 0.5834625 0.7253097 0.4579766 0.9950016 0.4701569
psi.1 <- AFMV.1$uniqueness
psi.2 <- AFMV.2$uniqueness
psi.1
## JP Citi Fargo Royal Exxon
## 0.4881494 0.2342357 0.5603556 0.8802602 0.9220490
psi.2
## JP Citi Fargo Royal Exxon
## 0.4165374 0.2746902 0.5420233 0.0050000 0.5298429
R.est.1 <- L.1%*%t(L.1)+diag(psi.1)
R.est.1
## 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
R.est.2 <- L.2%*%t(L.2)+diag(psi.2)
R.est.2
## 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
U.1 <- R - R.est.1
round(U.1,4)
## JP Citi Fargo Royal Exxon
## JP 0.0000 0.0062 0.0361 -0.1330 -0.0453
## Citi 0.0062 0.0000 -0.0061 0.0195 -0.0316
## Fargo 0.0361 -0.0061 0.0000 -0.0469 -0.0389
## Royal -0.1330 0.0195 -0.0469 0.0000 0.5868
## Exxon -0.0453 -0.0316 -0.0389 0.5868 0.0000
U.2 <- R - R.est.2
round(U.2,4)
## JP Citi Fargo Royal Exxon
## JP 0.0000 0.0000 -0.0026 -3e-04 0.0520
## Citi 0.0000 0.0000 0.0016 2e-04 -0.0331
## Fargo -0.0026 0.0016 0.0000 0e+00 0.0006
## Royal -0.0003 0.0002 0.0000 0e+00 0.0001
## Exxon 0.0520 -0.0331 0.0006 1e-04 0.0000
Conforme podemos observar a análise fatorial com dois fatores torna-se mais precisa do que com um, um resultado que já tinha sido observado acima, pois, a proporção da variabilidade explicada por dois fatores foi de: 64.6%.
ssc.1 <- sum2(L.1) # um fator
ssc.2 <- apply(L.2,2,"sum2") # dois fatores
ssc.1
## [1] 1.914953
ssc.2
## Factor1 Factor2
## 1.622061 1.609847
ssc.1/p # por um único fator
## [1] 0.3829907
ssc.2/p # por dois fatores
## Factor1 Factor2
## 0.3244121 0.3219693
prop.ac <- cumsum(ssc.2)/p # proporção da variabilidade acumulada m = 2
prop.ac
## Factor1 Factor2
## 0.3244121 0.6463815
# Testes de adequacao
nu.1 <- ((p-1)^2-p-1)/2
nu.2 <- ((p-2)^2-p-2)/2
nu.1;nu.2
## [1] 5
## [1] 1
m1 <- 1
W1 <- (n-1-(2*p + 4*m1 + 5)/6)*( log(det(R.est.1))-log(det(R)))
W1
## [1] 62.22057
Se o teste abaixo for verificado, podemos continuar que não é adequado a utilização da análise fatorial
W1 > qchisq(0.95,nu.1)
## [1] TRUE
valorp.1 <- 1-pchisq(W1,nu.1)
valorp.1
## [1] 4.221401e-12
Como a hipótese nula foi rejeitada, a quantidade de um fator não está adequada
m2 <- 2
W2 <- (n-1-(2*p + 4* m2 + 5)/6)*( log(det(R.est.2))-log(det(R)))
W2
## [1] 2.004672
W2 > qchisq(0.95,nu.2)
## [1] FALSE
valorp.2 <- 1-pchisq(W2,nu.2)
valorp.2
## [1] 0.1568152
Como a hipótese nula não foi rejeitada, a quantidade de dois fatores está adequada.
AFMV.2$PVAL
## objective
## 0.1600081
Comparado com o teste, vemos que dois fatores está bem adequado. Com um único fator, o resultado do p-valor foi de 4.22^{-12}. #### Estimação 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")
plot(AFMV.2.Reg$scores[,1],AFMV.2.Bar$scores[,2])
plot(AFMV.2.Reg$scores[,2],AFMV.2.Bar$scores[,1])
cor(AFMV.2.Reg$scores, AFMV.2.Bar$scores)
## Factor1 Factor2
## Factor1 1.000000e+00 1.840696e-17
## Factor2 3.429087e-17 1.000000e+00
AF <- prcomp(dados,scale=T)
delta <- AF$sdev^2
P <- AF$rotation
delta
## [1] 2.4372731 1.4070127 0.5005127 0.4000316 0.2551699
P
## PC1 PC2 PC3 PC4 PC5
## JP -0.4690832 0.3680070 -0.60431522 0.3630228 0.38412160
## Citi -0.5324055 0.2364624 -0.13610618 -0.6292079 -0.49618794
## Fargo -0.4651633 0.3151795 0.77182810 0.2889658 0.07116948
## Royal -0.3873459 -0.5850373 0.09336192 -0.3812515 0.59466408
## Exxon -0.3606821 -0.6058463 -0.10882629 0.4934145 -0.49755167
delta/p
## [1] 0.48745462 0.28140253 0.10010255 0.08000632 0.05103398
cumsum(delta)/p
## [1] 0.4874546 0.7688572 0.8689597 0.9489660 1.0000000
Ou seja, com dois fatores através de componentes principais temos: 76.89%.
ell.1 <- sqrt(delta[1])*P[,1]
ell.2 <- sqrt(delta[2])*P[,2]
L.1 <- ell.1 # um fator
L.2 <- cbind(ell.1,ell.2) # dois fatores
L.1
## JP Citi Fargo Royal Exxon
## -0.7323218 -0.8311791 -0.7262022 -0.6047155 -0.5630885
L.2
## ell.1 ell.2
## JP -0.7323218 0.4365209
## Citi -0.8311791 0.2804859
## Fargo -0.7262022 0.3738582
## Royal -0.6047155 -0.6939569
## Exxon -0.5630885 -0.7186401
Conforme já tinha sido encontrado através dos componentes principais estes dois escores, o primeiro uma soma ponderada e o segundo um contraste entre as ações dos bancos e as de petróleo.
h.1 <- L.1^2 # um fator
h.2 <- apply(L.2,1,"sum2") # dois fatores
h.1
## JP Citi Fargo Royal Exxon
## 0.5362953 0.6908587 0.5273697 0.3656808 0.3170686
h.2
## JP Citi Fargo Royal Exxon
## 0.7268458 0.7695311 0.6671396 0.8472571 0.8335122
psi_1 <- 1 - h.1
psi_2 <- 1 - h.2
psi_1
## JP Citi Fargo Royal Exxon
## 0.4637047 0.3091413 0.4726303 0.6343192 0.6829314
psi_2
## JP Citi Fargo Royal Exxon
## 0.2731542 0.2304689 0.3328604 0.1527429 0.1664878
R.est_1 <- L.1%*%t(L.1)+diag(psi_1)
R.est_2 <- L.2%*%t(L.2)+diag(psi_2)
R.est_1
## JP Citi Fargo Royal Exxon
## [1,] 1.0000000 0.6086906 0.5318137 0.4428464 0.4123620
## [2,] 0.6086906 1.0000000 0.6036041 0.5026269 0.4680274
## [3,] 0.5318137 0.6036041 1.0000000 0.4391457 0.4089161
## [4,] 0.4428464 0.5026269 0.4391457 1.0000000 0.3405083
## [5,] 0.4123620 0.4680274 0.4089161 0.3405083 1.0000000
R.est_2
## JP Citi Fargo Royal Exxon
## JP 1.00000000 0.7311286 0.6950107 0.1399197 0.09866056
## Citi 0.73112857 1.0000000 0.7084661 0.3079818 0.26645897
## Fargo 0.69501067 0.7084661 1.0000000 0.1797042 0.14024657
## Royal 0.13991966 0.3079818 0.1797042 1.0000000 0.83921362
## Exxon 0.09866056 0.2664590 0.1402466 0.8392136 1.00000000
U_1 <- R - R.est_2
round(U_1,3)
## JP Citi Fargo Royal Exxon
## JP 0.000 -0.099 -0.185 -0.025 0.056
## Citi -0.099 0.000 -0.134 0.014 -0.054
## Fargo -0.185 -0.134 0.000 0.003 0.006
## Royal -0.025 0.014 0.003 0.000 -0.156
## Exxon 0.056 -0.054 0.006 -0.156 0.000