library(readxl)
## Warning: package 'readxl' was built under R version 4.0.4
library(mvnTest)
## Warning: package 'mvnTest' was built under R version 4.0.5
## Loading required package: mvtnorm
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(MVN)
## Warning: package 'MVN' was built under R version 4.0.5
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
## sROC 0.1-2 loaded
data <- read_excel("D:/STIS/3SE3/12. APG/11. Paper UTS/Kelompok/Data_fix.xlsx")
data
## # A tibble: 82 x 5
## `Kota Inflasi` Normal Lebaran NormalCovid LebaranCovid
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 KOTA MEULABOH 136. 138. 106. 107.
## 2 KOTA BANDA ACEH 130. 130. 104. 105.
## 3 KOTA LHOKSEUMAWE 132. 132. 104. 104.
## 4 KOTA SIBOLGA 142. 146. 103. 103.
## 5 KOTA PEMATANG SIANTAR 141. 143. 103. 103.
## 6 KOTA MEDAN 143. 145. 103. 103.
## 7 KOTA PADANGSIDIMPUAN 134. 137. 105. 105.
## 8 KOTA PADANG 142. 143. 103. 104.
## 9 KOTA BUKITTINGGI 133. 134. 104. 104.
## 10 TEMBILAHAN 143. 144. 104. 105.
## # ... with 72 more rows
r1 = c(-1,-1,1,1)
r2 = c(-1,1,-1,1)
r3 = c(1,-1,-1,1)
C = as.matrix(rbind(r1,r2,r3))
y1 <- data$Normal
y2 <- data$Lebaran
y3 <- data$NormalCovid
y4 <- data$LebaranCovid
ybar1 <- mean(y1)
ybar2 <- mean(y2)
ybar3 <- mean(y3)
ybar4 <- mean(y4)
Ybar <- c(ybar1,ybar2,ybar3,ybar4)
Y <- matrix(c(y1,y2,y3,y4),nrow=length(y1),ncol=4)
Y
## [,1] [,2] [,3] [,4]
## [1,] 135.89 137.66 106.41 106.89
## [2,] 129.85 130.34 104.28 104.60
## [3,] 131.81 132.21 103.90 103.95
## [4,] 142.38 145.59 103.07 103.25
## [5,] 141.43 142.71 102.92 103.30
## [6,] 142.53 144.92 102.60 103.03
## [7,] 134.22 136.62 104.60 105.40
## [8,] 141.77 143.28 103.05 103.73
## [9,] 132.54 134.24 103.59 103.99
## [10,] 142.58 143.85 104.03 104.68
## [11,] 137.53 139.32 103.05 103.50
## [12,] 137.82 138.92 103.49 104.47
## [13,] 136.14 137.54 103.79 104.30
## [14,] 135.27 136.65 103.30 103.60
## [15,] 133.17 133.92 103.89 104.03
## [16,] 133.40 134.06 103.45 103.86
## [17,] 146.04 147.15 103.27 103.69
## [18,] 137.53 138.40 105.08 104.78
## [19,] 141.21 141.89 105.12 104.75
## [20,] 146.27 145.67 103.33 104.57
## [21,] 144.79 146.22 102.31 102.83
## [22,] 139.15 139.48 103.18 103.34
## [23,] 135.03 135.46 102.61 102.62
## [24,] 137.47 138.12 105.49 105.47
## [25,] 138.57 139.82 105.93 105.94
## [26,] 135.84 136.34 104.82 104.85
## [27,] 136.25 136.39 105.06 104.80
## [28,] 131.10 131.30 102.74 102.65
## [29,] 135.12 135.59 106.56 106.48
## [30,] 135.30 136.80 105.84 105.66
## [31,] 133.88 134.26 102.91 102.94
## [32,] 138.89 139.93 103.08 103.38
## [33,] 132.91 133.63 103.80 104.00
## [34,] 141.89 142.75 103.73 103.83
## [35,] 131.30 131.89 103.73 103.52
## [36,] 134.06 134.87 104.86 104.96
## [37,] 132.45 133.53 104.82 104.72
## [38,] 133.77 134.10 105.15 105.38
## [39,] 131.96 131.75 104.48 104.45
## [40,] 130.38 130.85 103.41 103.43
## [41,] 132.01 132.14 104.28 104.30
## [42,] 130.10 130.20 104.49 104.29
## [43,] 136.70 136.47 102.74 103.02
## [44,] 130.67 131.30 103.51 103.56
## [45,] 133.94 134.24 103.19 103.20
## [46,] 136.58 136.86 104.09 104.31
## [47,] 144.83 145.68 104.67 104.71
## [48,] 144.01 144.55 105.77 106.04
## [49,] 147.88 148.61 106.73 106.78
## [50,] 144.08 144.11 104.84 104.61
## [51,] 132.68 132.73 104.42 104.32
## [52,] 134.26 135.11 102.75 102.60
## [53,] 140.83 141.05 104.61 104.25
## [54,] 126.56 127.33 103.36 103.31
## [55,] 136.27 135.85 103.71 103.87
## [56,] 148.00 148.97 105.19 105.69
## [57,] 139.33 139.43 102.31 102.68
## [58,] 140.93 141.16 104.06 104.44
## [59,] 133.71 133.87 104.67 104.91
## [60,] 136.08 135.95 105.61 105.72
## [61,] 138.47 138.87 104.58 104.70
## [62,] 141.01 142.36 103.27 103.59
## [63,] 139.12 139.33 103.71 103.86
## [64,] 149.86 150.66 103.36 103.08
## [65,] 135.16 140.02 104.64 104.63
## [66,] 142.78 144.15 105.34 105.50
## [67,] 144.05 144.16 105.73 106.03
## [68,] 133.31 134.20 103.32 103.54
## [69,] 139.31 139.38 104.92 105.50
## [70,] 131.21 132.60 103.57 103.73
## [71,] 135.39 136.56 103.67 104.18
## [72,] 131.98 135.35 103.10 103.42
## [73,] 137.04 137.69 103.16 103.25
## [74,] 132.05 132.45 103.73 103.39
## [75,] 132.87 133.08 103.93 104.00
## [76,] 133.89 134.45 105.44 105.66
## [77,] 159.00 160.83 104.65 105.33
## [78,] 140.29 140.98 104.86 105.79
## [79,] 135.63 136.14 107.70 107.87
## [80,] 135.35 135.64 102.55 102.49
## [81,] 139.44 139.74 103.37 104.09
## [82,] 142.49 142.37 103.92 103.85
S <- cov(Y)
S
## [,1] [,2] [,3] [,4]
## [1,] 29.6351569 30.0728721 0.6076852 1.142606
## [2,] 30.0728721 31.2322569 0.5408383 1.126082
## [3,] 0.6076852 0.5408383 1.2268451 1.168823
## [4,] 1.1426061 1.1260824 1.1688228 1.210026
n <- length(y1)
q <- ncol(Y)
T2 <- n*t(C%*%Ybar)%*%solve(C%*%S%*%t(C))%*%(C%*%Ybar)
T2
## [,1]
## [1,] 3305.324
Ftable = qf(1-0.05, q-1,n-q+1)
C2 = Ftable*(n-1)*(q-1)/(n-q+1)
C2
## [1] 8.367396
CY = C%*%Ybar
CSC = C%*%S%*%t(C)
{
a <- diag(3)
p = 3
LowerBound <- c()
UpperBound <- c()
for (i in 1:p) {
LowerBound[i] <- (t(a[,i])%*%CY) - sqrt(C2*t(a[,i])%*%CSC%*%a[,i]/n)
UpperBound[i] <- (t(a[,i])%*%CY) + sqrt(C2*t(a[,i])%*%CSC%*%a[,i]/n)
}
cat("95% simultaneous confidence interval for each variable:\n\n",
LowerBound[1],"< mu1 <",UpperBound[1],"\n",
LowerBound[2],"< mu2 <",UpperBound[2],"\n",
LowerBound[3],"< mu3 <",UpperBound[3])
}
## 95% simultaneous confidence interval for each variable:
##
## -70.69521 < mu1 < -63.72723
## 0.6867627 < mu2 < 1.300067
## -0.8874413 < mu3 < -0.3452416