###################################################
### chunk number 1: 
###################################################
options(show.signif.stars=FALSE,digits=4,width=80)


###################################################
### chunk number 2: 
###################################################
sxx <- matrix(c(1.,.7,.7,1),nrow=2)
syy <- matrix(c(1,.6,.6,1),nrow=2)
sxy <- matrix(c(.5,.6,.3,.4),nrow=2,byrow=T)
(R <- cbind (rbind(sxx,sxy), rbind(t(sxy),syy)))
eigen(R)$values


###################################################
### chunk number 3: 
###################################################
eigen.sxx <- eigen(sxx)
eigen.syy <- eigen(syy)
(sxxinvsqrt <- eigen.sxx$vectors %*% diag(1/sqrt(eigen.sxx$values)) %*%
              t(eigen.sxx$vectors))
(syyinvsqrt <- eigen.syy$vectors %*% diag(1/sqrt(eigen.syy$values)) %*%
              t(eigen.syy$vectors))


###################################################
### chunk number 4: 
###################################################
(mat <- sxxinvsqrt %*% sxy %*% syyinvsqrt)
(sv <- svd(mat))


###################################################
### chunk number 5: 
###################################################
t(t(sv$u) %*% sxxinvsqrt)
t(t(sv$v) %*% syyinvsqrt)


###################################################
### chunk number 6: one
###################################################
pairs(LifeCycleSavings[,c(2,3,1,4,5)])


###################################################
### chunk number 7: 
###################################################
pop <- LifeCycleSavings[, 2:3]
oec <- LifeCycleSavings[, -(2:3)]
(c1<-cancor(pop, oec))


###################################################
### chunk number 8: 
###################################################
(c2<- cancor(scale(pop),scale(oec)))


###################################################
### chunk number 9: two
###################################################
plot(as.matrix(pop) %*% c1$xcoef[,1],
     as.matrix(scale(pop))%*% c2$xcoef[,1],xlab="Raw",ylab="Standardized")


###################################################
### chunk number 10: three
###################################################
d <- cbind(as.matrix(pop) %*% c1$xcoef,as.matrix(oec) %*% c1$ycoef[,1:2])
colnames(d) <- c("X1", "X2", "Y1", "Y2")
pairs(d)


###################################################
### chunk number 11: four
###################################################
plot( as.matrix(pop) %*% c1$xcoef[,1],
      as.matrix(oec) %*% c1$ycoef[,1],type="n",
      xlab="X-canonical vector",ylab="Y-canonical vector")
text( as.matrix(pop) %*% c1$xcoef[,1],
      as.matrix(oec) %*% c1$ycoef[,1],substr(rownames(oec),1,8),cex=.6)


###################################################
### chunk number 12: 
###################################################
n <- dim(pop)[1]
p <- 2
q <- 3
(test <- -(n-(p+q+3)/2) *log(prod(1-c1$cor^2)))
pchisq(test,p*q,lower.tail=FALSE)


###################################################
### chunk number 13: 
###################################################
(test <- -(n-(p+q+3)/2) *log(prod(1-c1$cor[-1]^2)))
pchisq(test,(p-1)*(q-1),lower.tail=FALSE)


###################################################
### chunk number 14: 
###################################################
library("systemfit")
data(Kmenta)
eqDemand <- consump ~ price + income
eqSupply <- consump ~ price + farmPrice + trend
eqSystem <- list(demand = eqDemand, supply = eqSupply)
fitols <- systemfit(eqSystem, data=Kmenta)
summary(fitsur <- systemfit(eqSystem, method = "SUR", data=Kmenta))
compareCoefs(fitols, fitsur)