Purpose - Johnson Inference of Means
Again, consider the monthly log returns of the CRSP ten decile indices from January 2000 to December 2009. Construct

  1. Benferroni, and
  2. marginal confidence intervals for the means of the ten portfolios.
  3. compute the lengths of confidence intervals.
  4. simultaneous T2 intervals
> folder <- "C:/Cauldron/Books/Statistics/Multivariate/TSAY/"
> x <- read.table(paste(folder, "m-dec0009.txt", sep = ""), header = T)
> x <- x[, -1]
> x <- log(1 + x)
> mu <- colMeans(x)
> S <- cov(x)
> n <- dim(x)[1]
> p <- ncol(x)
> S.ind <- diag(S)/n
> alpha <- 0.05
> interval <- sqrt(((n - 1) * p * qf(1 - alpha, p, n - p)/(n -
+     p)) * (S.ind))
> cbind(mu - interval, mu + interval)
             [,1]       [,2]
cap1  -0.02025303 0.04990609
cap2  -0.01908871 0.03849229
cap3  -0.01935686 0.03371115
cap4  -0.01843208 0.03485084
cap5  -0.02066861 0.03326123
cap6  -0.02233284 0.03209764
cap7  -0.02246355 0.03272190
cap8  -0.02264752 0.02938323
cap9  -0.02063713 0.02781794
cap10 -0.02096751 0.01983116
  1. Bonferroni Confidence Intervals
> interval <- qt(1 - (alpha/(2 * p)), n - 1) * sqrt(S.ind)
> cbind(mu - interval, mu + interval)
              [,1]       [,2]
cap1  -0.007202043 0.03685510
cap2  -0.008377507 0.02778108
cap3  -0.009485158 0.02383945
cap4  -0.008520400 0.02493916
cap5  -0.010636589 0.02322921
cap6  -0.012207696 0.02197249
cap7  -0.012197959 0.02245631
cap8  -0.012968774 0.01970448
cap9  -0.011623526 0.01880434
cap10 -0.013378153 0.01224180
  1. Individial Confidence Intervals
> interval <- qt(1 - (alpha/(2)), n - 1) * sqrt(S.ind)
> cbind(mu - interval, mu + interval)
               [,1]       [,2]
cap1  -0.0004231019 0.03007616
cap2  -0.0028138930 0.02221747
cap3  -0.0043576007 0.01871189
cap4  -0.0033720771 0.01979084
cap5  -0.0054257600 0.01801838
cap6  -0.0069484926 0.01671329
cap7  -0.0068658100 0.01712416
cap8  -0.0079414385 0.01467714
cap9  -0.0069416820 0.01412249
cap10 -0.0094360878 0.00829973

4.Asymptotic Confidence

> interval <- sqrt(qchisq((1 - alpha), p) * S.ind)
> cbind(mu - interval, mu + interval)
             [,1]       [,2]
cap1  -0.01812543 0.04777850
cap2  -0.01734255 0.03674613
cap3  -0.01774756 0.03210185
cap4  -0.01681626 0.03323502
cap5  -0.01903317 0.03162579
cap6  -0.02068222 0.03044702
cap7  -0.02079003 0.03104838
cap8  -0.02106968 0.02780538
cap9  -0.01916771 0.02634852
cap10 -0.01973028 0.01859393

The same intervals can be generated by the cregion code from Johnson.

> source(paste(folder, "cregion.R", sep = ""))
> confreg(x)
[1] "C.R. based on T^2"
             [,1]       [,2]
 [1,] -0.02025303 0.04990609
 [2,] -0.01908871 0.03849229
 [3,] -0.01935686 0.03371115
 [4,] -0.01843208 0.03485084
 [5,] -0.02066861 0.03326123
 [6,] -0.02233284 0.03209764
 [7,] -0.02246355 0.03272190
 [8,] -0.02264752 0.02938323
 [9,] -0.02063713 0.02781794
[10,] -0.02096751 0.01983116
[1] "CR based on individual t"
               [,1]       [,2]
 [1,] -0.0004231019 0.03007616
 [2,] -0.0028138930 0.02221747
 [3,] -0.0043576007 0.01871189
 [4,] -0.0033720771 0.01979084
 [5,] -0.0054257600 0.01801838
 [6,] -0.0069484926 0.01671329
 [7,] -0.0068658100 0.01712416
 [8,] -0.0079414385 0.01467714
 [9,] -0.0069416820 0.01412249
[10,] -0.0094360878 0.00829973
[1] "CR based on Bonferroni"
              [,1]       [,2]
 [1,] -0.007202043 0.03685510
 [2,] -0.008377507 0.02778108
 [3,] -0.009485158 0.02383945
 [4,] -0.008520400 0.02493916
 [5,] -0.010636589 0.02322921
 [6,] -0.012207696 0.02197249
 [7,] -0.012197959 0.02245631
 [8,] -0.012968774 0.01970448
 [9,] -0.011623526 0.01880434
[10,] -0.013378153 0.01224180
[1] "Asymp. simu. CR"
             [,1]       [,2]
 [1,] -0.01812543 0.04777850
 [2,] -0.01734255 0.03674613
 [3,] -0.01774756 0.03210185
 [4,] -0.01681626 0.03323502
 [5,] -0.01903317 0.03162579
 [6,] -0.02068222 0.03044702
 [7,] -0.02079003 0.03104838
 [8,] -0.02106968 0.02780538
 [9,] -0.01916771 0.02634852
[10,] -0.01973028 0.01859393

Basically once you code you understand that calculations as simple as getting a confidence interval for a p dim vector is complicated.