## Question 8.10

## part(a)

setwd("/Users/pszhong/Documents/MSU-Teaching/STT-843-Spring-2018")
stockdata<-read.table(file="stocks-dat.txt")
colnames(stockdata)<-c("JPMorgan","Citibank","WellsFargo","RoyalDutchShell","ExxonMobil")
Sn<-cov(stockdata)
emat<-eigen(Sn)
ehats<-emat$vectors

## part(b)

evals<-emat$values
prop3<-sum(evals[1:3])/sum(evals)

## part (c)

n<-dim(stockdata)[1]
cbind(evals-qnorm(1-0.1/6)*sqrt(2*evals^2/n),evals+qnorm(1-0.1/6)*sqrt(2*evals^2/n))

## Question 8.28

## part (a)

setwd("/Users/pszhong/Documents/MSU-Teaching/STT-843-Spring-2018")
malidata<-read.table(file="mali-dat.txt",header=TRUE)
plot(malidata$Family,malidata$DistRD,xlab="Family", ylab="DistRD")
rowoutlier<-c(1:nrow(malidata))[malidata$Family>140]
coloutlier<-c(1:nrow(malidata))[malidata$DistRD>400]
outlier<-c(rowoutlier,coloutlier)
text(malidata$Family[outlier]+0.5, malidata$DistRD[outlier]-14, labels=outlier,cex=0.7)

plot(malidata$DistRD,malidata$Cattle,xlab="DistRD",ylab="Cattle")
rowoutlier1<-c(1:nrow(malidata))[malidata$DistRD>400]
coloutlier1<-c(1:nrow(malidata))[malidata$Cattle>100]
outlier1<-c(rowoutlier1,coloutlier1)
text(malidata$DistRD[coloutlier1]+0.5, malidata$Cattle[coloutlier1]-3, labels=coloutlier1,cex=0.7)
text(malidata$DistRD[rowoutlier1[1]]+10, malidata$Cattle[rowoutlier1[1]]+3,labels=rowoutlier1[1],cex=0.7)
text(malidata$DistRD[rowoutlier1[2]]-10, malidata$Cattle[rowoutlier1[2]]+3,labels=rowoutlier1[2],cex=0.7)

malidataNoOutliers<-malidata[-unique(c(outlier,outlier1)),]

## part (b)

mali.pca<-prcomp(malidataNoOutliers,center=TRUE,scale=TRUE)
## Plot the variances of the principal components
plot(mali.pca,type="l",main="Scree plot")
drops<-diff((mali.pca$sdev)^2)
## Variance explained by each PCs
totalvariance<-sum((mali.pca$sdev)^2)
varianceratio<-(mali.pca$sdev)^2/totalvariance
plot(varianceratio,xlab="Component number")


## Question 10.2

## part (a)

Sigma11<-matrix(c(8,2,2,5),2,2)
Sigma12<-matrix(c(3,-1,1,3),2,2)
Sigma21<-matrix(c(3,1,-1,3),2,2)
Sigma22<-matrix(c(6,-2,-2,7),2,2)
C0<-solve(Sigma11)%*%Sigma12%*%solve(Sigma22)%*%Sigma21
eigen(C0)

## part (b)

InvSigmaHalf<-function(Sigma)
{
EigenSigma<-eigen(Sigma)
InvSighalf<-EigenSigma$vectors%*%diag(1/sqrt(EigenSigma$values))%*%t(EigenSigma$vectors)
return(InvSighalf)
}
C1<-InvSigmaHalf(Sigma11)%*%Sigma12%*%solve(Sigma22)%*%Sigma21%*%InvSigmaHalf(Sigma11)
EC1<-eigen(C1)
C2<-InvSigmaHalf(Sigma22)%*%Sigma21%*%solve(Sigma11)%*%Sigma12%*%InvSigmaHalf(Sigma22)
EC2<-eigen(C2)
U12<-t(EC1$vectors)%*%InvSigmaHalf(Sigma11)
V12<-t(EC2$vectors)%*%InvSigmaHalf(Sigma22)

## part (c)

mu<-c(-3,2,0,1)
A<-rbind(cbind(U12,matrix(0,2,2)),cbind(matrix(0,2,2),V12))
Emus<-A%*%mu
Sigma<-rbind(cbind(Sigma11,Sigma12),cbind(Sigma21,Sigma22))
CovUVs<-A%*%Sigma%*%t(A)
round(CovUVs,3)

## Question 10.10

## part (a)

R11<-matrix(c(1,0.615,0.615,1),2,2)
R12<-matrix(c(-0.111,-0.195,-0.266,-0.085),2,2)
R21<-matrix(c(-0.111,-0.266,-0.195,-0.085),2,2)
R22<-matrix(c(1,-0.269,-0.269,1),2,2)
RR<-solve(R11)%*%R12%*%solve(R22)%*%R21
eigen(RR)

## part (b)

RR1<-InvSigmaHalf(R11)%*%R12%*%solve(R22)%*%R21%*%InvSigmaHalf(R11)
ERs<-eigen(RR1)
fu1<-t(ERs$vectors)%*%InvSigmaHalf(R11)
RR2<-InvSigmaHalf(R22)%*%R21%*%solve(R11)%*%R12%*%InvSigmaHalf(R22)
ERs2<-eigen(RR2)
fu2<-t(ERs2$vectors)%*%InvSigmaHalf(R22)

## Question 10.13

## part (a)

R11<-diag(rep(1,5))
R11[lower.tri(R11)]<-c(.754,-.690,-0.446,0.692,-.712,-.515,.412,0.323,-.444,-0.334)
R11<-R11+t(R11)-diag(rep(1,5))
R21<-matrix(c(-0.605,-.479,.780,-.152,-.722,-0.419,0.542,-.102,0.737,0.361,-0.546,0.172,.527,0.461,-0.393,-.019,-.383,-.505,.737,-0.148),4,5)
R22<-diag(rep(1,4))
R22[lower.tri(R22)]<-c(.251,-.490,0.250,-.434,-.079,-.163)
R22<-R22+t(R22)-diag(rep(1,4))
RR0<-solve(R22)%*%R21%*%solve(R11)%*%t(R21)
rhostars2<-eigen(RR0)
n<-138
p<-5
q<-4
Sigma12teststat<- -(n-1-(p+q+1)/2)*log(prod(1-rhostars2$values))
Sigma12testpval<-1-pchisq(Sigma12teststat,p*q)
rho1teststat<- -(n-1-(p+q+1)/2)*log(prod(1-rhostars2$values[2:4]))
rho1testpval<-1-pchisq(rho1teststat,(p-1)*(q-1))
rho2teststat<- -(n-1-(p+q+1)/2)*log(prod(1-rhostars2$values[3:4]))
rho2testpval<-1-pchisq(rho2teststat,(p-2)*(q-2))

## part (b)

RR1<-InvSigmaHalf(R11)%*%t(R21)%*%solve(R22)%*%R21%*%InvSigmaHalf(R11)
Us<-t(eigen(RR1)$vectors)%*%InvSigmaHalf(R11)
RR2<-InvSigmaHalf(R22)%*%R21%*%solve(R11)%*%t(R21)%*%InvSigmaHalf(R22)
Vs<-t(eigen(RR2)$vectors)%*%InvSigmaHalf(R22)

## part (c)

invU<-solve(Us)
U1prop<-sum(invU[,1]^2)/p
invV<-solve(Vs)
V1prop<-sum(invV[,1]^2)/q