R : Copyright 2003, The R Development Core Team Version 1.6.2 (2003-01-10) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type `license()' or `licence()' for distribution details. R is a collaborative project with many contributors. Type `contributors()' for more information. Type `demo()' for some demos, `help()' for on-line help, or `help.start()' for a HTML browser interface to help. Type `q()' to quit R. > invisible(options(echo = TRUE)) > library(cluster) Loading required package: mva > > rally <- read.table("rally.csv",header=T,sep=",") > source("psopts.r") > > n <- ncol(rally) > X <- t(rally[,2:n]) > n <- nrow(X) > p <- ncol(X) > > X.pca <- princomp(X,covmat = cov.wt(X,center=F)) > #X.pca <- princomp(X) > print(summary(X.pca)) Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Standard deviation 16.2640262 4.78800399 2.65696871 2.37510859 2.19239962 Proportion of Variance 0.8270687 0.07167942 0.02207285 0.01763812 0.01502882 Cumulative Proportion 0.8270687 0.89874811 0.92082095 0.93845908 0.95348790 Comp.6 Comp.7 Comp.8 Comp.9 Standard deviation 2.02938671 1.84673916 1.619369616 1.357286279 Proportion of Variance 0.01287701 0.01066342 0.008199312 0.005760078 Cumulative Proportion 0.96636491 0.97702833 0.985227642 0.990987720 Comp.10 Comp.11 Standard deviation 1.274572973 1.121530431 Proportion of Variance 0.005079429 0.003932852 Cumulative Proportion 0.996067148 1.000000000 > W <- -X.pca$loadings > X.px <- X%*%W > > postscript(file="rally_var.eps") > screeplot(X.pca,main="") > dev.off() null device 1 > > postscript(file="rally_pca.eps") > plot(X.px[,1:2],type="n") > text(X.px[,1:2],labels = rownames(X.px),cex=0.5) > dev.off() null device 1 > > list = c("ward","complete","average","single") > > for (method in list) { + h <- hclust(dist(X),method=method) + print(summary(h)) + print(cutree(h,k=3)) + meth <- substr(method,1,4) + filename <- paste("rally_hcl_",meth,".eps",sep="") + postscript(file=filename) + plot(h,cex=0.6,lwd=0.6,ann=F) + dev.off() + } Length Class Mode merge 156 -none- numeric height 78 -none- numeric order 79 -none- numeric labels 79 -none- character method 1 -none- character call 3 -none- call dist.method 1 -none- character sp500 nyse nasdaq S01 S02 S10 S12 S13 S14 S15 S16 1 1 2 3 3 1 1 1 1 1 1 S17 S20 S21 S22 S23 S24 S25 S26 S27 S28 S29 1 1 1 1 1 1 1 1 1 1 1 S30 S31 S32 S33 S34 S35 S36 S37 S38 S39 S40 1 1 2 1 1 2 2 1 1 1 1 S41 S42 S44 S45 S47 S48 S49 S50 S51 S52 S53 2 1 1 1 1 1 1 1 1 2 1 S54 S55 S56 S57 S58 S59 S60 S61 S62 S63 S64 1 1 2 2 1 1 1 1 2 1 1 S65 S67 S70 S72 S73 S75 S78 S79 S80 S82 S89 1 1 1 1 2 1 2 1 1 1 1 tcm10y tcm1y tcm3y tcm5y tcm7y yen euro SL SM SH BL 3 3 3 3 3 3 3 1 1 1 1 BM BH 1 1 Length Class Mode merge 156 -none- numeric height 78 -none- numeric order 79 -none- numeric labels 79 -none- character method 1 -none- character call 3 -none- call dist.method 1 -none- character sp500 nyse nasdaq S01 S02 S10 S12 S13 S14 S15 S16 1 1 2 3 3 1 1 1 1 1 1 S17 S20 S21 S22 S23 S24 S25 S26 S27 S28 S29 1 1 1 1 2 1 1 1 1 1 1 S30 S31 S32 S33 S34 S35 S36 S37 S38 S39 S40 1 1 2 1 1 2 2 1 1 1 1 S41 S42 S44 S45 S47 S48 S49 S50 S51 S52 S53 2 1 1 1 1 1 1 1 1 2 2 S54 S55 S56 S57 S58 S59 S60 S61 S62 S63 S64 1 2 2 2 1 1 1 1 2 1 1 S65 S67 S70 S72 S73 S75 S78 S79 S80 S82 S89 1 1 2 1 2 1 2 2 1 1 1 tcm10y tcm1y tcm3y tcm5y tcm7y yen euro SL SM SH BL 3 3 3 3 3 3 3 1 1 1 1 BM BH 1 1 Length Class Mode merge 156 -none- numeric height 78 -none- numeric order 79 -none- numeric labels 79 -none- character method 1 -none- character call 3 -none- call dist.method 1 -none- character sp500 nyse nasdaq S01 S02 S10 S12 S13 S14 S15 S16 1 1 2 3 3 1 1 1 1 1 1 S17 S20 S21 S22 S23 S24 S25 S26 S27 S28 S29 1 1 1 1 1 1 1 1 1 1 1 S30 S31 S32 S33 S34 S35 S36 S37 S38 S39 S40 1 1 2 1 1 2 2 1 1 1 1 S41 S42 S44 S45 S47 S48 S49 S50 S51 S52 S53 2 1 1 1 1 2 1 2 1 2 1 S54 S55 S56 S57 S58 S59 S60 S61 S62 S63 S64 1 1 2 2 1 1 1 1 2 1 1 S65 S67 S70 S72 S73 S75 S78 S79 S80 S82 S89 1 1 1 2 2 1 2 1 1 1 1 tcm10y tcm1y tcm3y tcm5y tcm7y yen euro SL SM SH BL 3 3 3 3 3 3 3 1 1 1 1 BM BH 1 1 Length Class Mode merge 156 -none- numeric height 78 -none- numeric order 79 -none- numeric labels 79 -none- character method 1 -none- character call 3 -none- call dist.method 1 -none- character sp500 nyse nasdaq S01 S02 S10 S12 S13 S14 S15 S16 1 1 1 1 2 1 1 1 1 1 1 S17 S20 S21 S22 S23 S24 S25 S26 S27 S28 S29 1 1 1 1 1 1 1 1 1 1 1 S30 S31 S32 S33 S34 S35 S36 S37 S38 S39 S40 1 1 1 1 1 1 1 1 1 1 1 S41 S42 S44 S45 S47 S48 S49 S50 S51 S52 S53 3 1 1 1 1 1 1 1 1 1 1 S54 S55 S56 S57 S58 S59 S60 S61 S62 S63 S64 1 1 1 1 1 1 1 1 1 1 1 S65 S67 S70 S72 S73 S75 S78 S79 S80 S82 S89 1 1 1 1 1 1 1 1 1 1 1 tcm10y tcm1y tcm3y tcm5y tcm7y yen euro SL SM SH BL 1 1 1 1 1 1 1 1 1 1 1 BM BH 1 1 > > rn <- rownames(X) # Check that column numbers > print(c(rn[70],rn[1],rn[3])) # are correct. [1] "tcm5y" "sp500" "nasdaq" > > tcm5y <- X[70,] > sp500 <- X[1,] > nasdaq <- X[3,] > > rn.px <- rownames(X.px) # Check that column numbers > print(c(rn.px[70],rn.px[1],rn.px[3])) # are correct. [1] "tcm5y" "sp500" "nasdaq" > > tcm5y.px <- X.px[70,1:2] > sp500.px <- X.px[1,1:2] > nasdaq.px <- X.px[3,1:2] > > ctr <- rbind(tcm5y, sp500, nasdaq) > ctr.px <- ctr%*%W > > postscript(file="rally_km0.eps") > plot(X.px[,1:2], type="n") > text(X.px[,1:2],labels = rownames(X),cex=0.5) > points(ctr.px[,1:2],pch=19,cex=1.5,col=c("black","red","green")) > dev.off() null device 1 > > clust <- mat.or.vec(n,3) > for (i in 1:n) { + clust[i,] <- + order(c(sum((X[i,]-tcm5y)^2), sum((X[i,]-sp500)^2), sum((X[i,]-nasdaq)^2))) + } > > postscript(file="rally_km1.eps") > plot(X.px[,1:2], type="n") > text(X.px[,1:2], labels = rownames(X),cex=0.5,col=clust[,1]) > points(ctr.px[,1:2],pch=19,cex=1.5,col=c("black","red","green")) > dev.off() null device 1 > > ctr <- rbind(tcm5y, sp500, nasdaq) > ctr.px <- ctr%*%W > > km <- kmeans(X,ctr,20) > print(summary(km)) Length Class Mode cluster 79 -none- numeric centers 33 -none- numeric withinss 3 -none- numeric size 3 -none- numeric > > ctr <- km$centers > ctr.px <- ctr%*%W > > postscript(file="rally_km2.eps") > plot(X.px[,1:2], type="n") > text(X.px[,1:2], labels = rownames(X),cex=0.5,col=km$cluster) > points(ctr.px[,1:2],pch=19,cex=1.5,col=c("black","red","green")) > dev.off() null device 1 > > postscript(file="rally_km.eps") > plot(X.px[,1:2], type="n") > text(X.px[,1:2], labels = rownames(X),cex=0.5,col=km$cluster) > dev.off() null device 1 > > postscript(file="rally_ret.eps") > plot(x=X[,3],y=X[,9], type="n", + xlab="October 21, 1987, returns", + ylab="July, 24, 2002, returns" + ) > text(x=X[,3],y=X[,9], labels = rownames(X),cex=0.5,col=km$cluster) > dev.off() null device 1 > proc.time() [1] 0.51 0.05 0.82 0.00 0.00 >