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 > > crash <- read.table("crash.csv",header=T,sep=",") > source("psopts.r") > > n <- ncol(crash) > X <- t(crash[,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 23.6338836 4.05770733 2.69276802 2.248331709 1.861033732 Proportion of Variance 0.9342069 0.02753812 0.01212749 0.008454611 0.005792704 Cumulative Proportion 0.9342069 0.96174501 0.97387250 0.982327108 0.988119812 Comp.6 Comp.7 Comp.8 Comp.9 Standard deviation 1.650114015 1.432621693 1.000160540 0.894108329 Proportion of Variance 0.004554081 0.003432701 0.001673063 0.001337067 Cumulative Proportion 0.992673894 0.996106594 0.997779657 0.999116724 Comp.10 Standard deviation 0.7267109619 Proportion of Variance 0.0008832757 Cumulative Proportion 1.0000000000 > W <- -X.pca$loadings > X.px <- X%*%W > > postscript(file="crash_var.eps") > screeplot(X.pca,main="") > dev.off() null device 1 > > postscript(file="crash_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("crash_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 1 2 2 2 2 1 2 1 1 S17 S20 S21 S22 S23 S24 S25 S26 S27 S28 S29 1 1 2 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 1 2 2 1 1 1 2 1 2 1 1 S54 S55 S56 S57 S58 S59 S60 S61 S62 S63 S64 1 1 1 1 2 1 1 1 1 2 2 S65 S67 S70 S72 S73 S75 S78 S79 S80 S82 S89 2 2 1 1 1 1 1 1 1 2 2 tcm10y tcm1y tcm3y tcm5y tcm7y yen euro SL SM SH BL 3 3 3 3 3 3 3 1 2 2 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 2 2 2 2 1 2 1 1 S17 S20 S21 S22 S23 S24 S25 S26 S27 S28 S29 1 1 2 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 1 2 2 1 1 1 2 1 2 1 1 S54 S55 S56 S57 S58 S59 S60 S61 S62 S63 S64 1 1 1 1 2 1 1 1 1 2 2 S65 S67 S70 S72 S73 S75 S78 S79 S80 S82 S89 2 2 1 1 1 1 1 1 1 2 2 tcm10y tcm1y tcm3y tcm5y tcm7y yen euro SL SM SH BL 3 3 3 3 3 3 3 1 2 2 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 1 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 2 1 1 1 2 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 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 1 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 3 3 3 3 3 3 3 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="crash_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="crash_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 30 -none- numeric withinss 3 -none- numeric size 3 -none- numeric > > ctr <- km$centers > ctr.px <- ctr%*%W > > postscript(file="crash_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="crash_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="crash_ret.eps") > plot(x=X[,3],y=X[,8], type="n", + xlab="October 19, 1987, returns", + ylab="August 31, 1998, 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.53 0.02 0.56 0.00 0.00 >