Purpose - Johnson Inference of Means

  1. Consider the monthly simple returns of CRSP ten decile indices from January 2000 to December 2009. The portfolios are rebalanced annually and returns are obtained from CRSP. Transform the simple returns into log returns. Answer the following questions:

Obtain scatterplots of pairs of the decile log returns. Comments on the relation-
ship between the monthly log returns.

> folder <- "C:/Cauldron/Books/Statistics/Multivariate/TSAY/"
> x <- read.table(paste(folder, "m-dec0009.txt", sep = ""), header = T)
> x <- x[, -1]
> x.ln <- log(1 + x)
> plot(x)

DP_2-001.jpg

Are the monthly log returns jointly normal? You may use QQ-plot.

> par(mfrow = c(3, 2))
> for (i in 1:6) {
+     qqnorm(x.ln[, i], main = i)
+     qqline(x.ln[, i])
+ }

DP_2-002.jpg

> par(mfrow = c(2, 2))
> for (i in 7:10) {
+     qqnorm(x.ln[, i], main = i)
+     qqline(x.ln[, i])
+ }

DP_2-003.jpg

Looks like they satisfy marginal normality

To test whether they are jointly normal or not

> library(mvoutlier)
> par(mfrow = c(1, 1))
> chisq.plot(x.ln[, 1:10], ask = F)
$outliers
NULL
  1. Correlation of the chi square plot
> qqchi2(x)
[1] "correlation coefficient:"
[1] 0.9496503

DP_2-005.jpg

Assume that the monthly log returns are jointly normal. Test Ho mu= 0 versus Ha : m not eq to 0, Perform the test and draw your conclusion.

> mu <- colMeans(x)
> n <- dim(x)[1]
> p <- ncol(x)
> S <- cov(x)
> S.inv <- solve(S)
> alpha <- 0.1
> Hotelling <- n * t(mu) %*% S.inv %*% mu
> critical.val <- (n - 1) * p * qf(1 - alpha, p, n - p)/(n - p)

Cannot reject the null.
Ok, one thing I dont understand is that the chisquare plot from qqchi2 is different mvoutlier package. Want to investigate it closely

> da <- x
> nr <- dim(da)[1]
> nc <- dim(da)[2]
> cm <- matrix(colMeans(da), 1, nc)
> o <- matrix(rep(1, nr), nr, 1)
> dev <- da - kronecker(o, cm)
> dev <- as.matrix(dev)
> s <- t(dev) %*% dev/(nr - 1)
> si <- solve(s)
> d2 <- sort(diag(dev %*% si %*% t(dev)))
> prob <- (c(1:nr) - 0.5)/nr
> q1 <- qchisq(prob, nc)
> plot(q1, d2, xlab = "Quantile of chi-square", ylab = "d")
> fit <- lsfit(q1, d2)
> fitted <- d2 - fit$residuals
> lines(q1, fitted)

DP_2-007.jpg

chisq.plot in outliers package

> quan <- 1/2
> covr <- covMcd(x, alpha = quan)
> dist <- mahalanobis(x, center = covr$center, cov = covr$cov)
> s <- sort(dist, index = TRUE)
> q <- (0.5:length(dist))/length(dist)
> qchi <- qchisq(q, df = ncol(x))
> plot(s$x, qchi, xlab = "Ordered robust MD^2", ylab = "Quantiles of Chi_p^2",
+     main = "Chi^2-Plot", col = 3)

DP_2-008.jpg

I think I wll go with chisquare plot to test the joint normality and use chisq to look at outliers .