Redo Midterm Exam

Name : Violetta Ammara Irianti Agung

Student ID : 5003251126

Question 1 – Dataset: AirPassengers

df <- matrix(AirPassengers,nrow=12,ncol=12)
colnames(df) <- 1949:1960
rownames(df) <- month.abb

# (a) arith_mean manual
arith_mean <- function(x){
  n <- 0
  total <- 0
  for(val in x){
    total <- total+val
    n <- n+1
  }
  return(total/n)
}

# apply arith_mean() to each column of df
Y_bar <- numeric(12)
for (j in 1:12){
  Y_bar[j] <- arith_mean(df[,j])
}

# display results
names(Y_bar) <- 1949:1960
print(Y_bar)

    1949     1950     1951     1952     1953     1954 
126.6667 139.6667 170.1667 197.0000 225.0000 238.9167 
    1955     1956     1957     1958     1959     1960 
284.0000 328.2500 368.4167 381.0000 428.3333 476.1667 

# validation
print(colMeans(df)==Y_bar)

1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE 

# (b) geo_mean_growth()
geo_mean_growth <- function(y){
  n <- 0
  for(val in y) n <- n+1

  # (i) growth factors r_t = Y_t / Y_{t-1}, t = 2..n
  r <- numeric(n-1)
  for(t in 2:n){
    r[t-1] <- y[t]/y[t-1]
  }

  # (ii) geometric mean GM = (r2 * r3 * ... * rn)^(1/(n-1))
  product <- 1
  for(val in r){
    product <- product*val
  }
  GM <- product^(1/(n-1))

  # (iii) convert to percentage
  g_bar <- (GM-1)*100

  return(list(r=r,
              GM=GM,
              g_bar=g_bar))
}

# apply to Y_bar
results_GM <- geo_mean_growth(Y_bar)
cat(sprintf("Geometric mean growth rate: %.4f%%\n", results_GM$g_bar))

Geometric mean growth rate: 12.7928%

# (c) arithmetic mean of growth factors r_t

# arithmetic mean of r_t
AM <- arith_mean(results_GM$r)

# convert into percentage growth
g_bar_AM <- (AM-1)*100
cat(sprintf("Arithmetic mean growth rate: %.4f%%\n", g_bar_AM))

Arithmetic mean growth rate: 12.9056%

# (d) predict Y_bar_1960 and compare with the actual value

# extract from Y_bar
Y_1949 <- Y_bar["1949"]
Y_1960 <- Y_bar["1960"]

# prediction using geometric mean growth rate
Y_hat_GM <- Y_1949*(1+(results_GM$g_bar/100))^11

# prediction using arithmetic mean growth rate
Y_hat_AM <- Y_1949*(1+(g_bar_AM/100))^11

# absolute prediction error for each |Y_hat_1960 - Y_bar_1960|
err_GM <- abs(Y_hat_GM-Y_1960)
err_AM <- abs(Y_hat_AM-Y_1960)

{
  cat(sprintf("Y_hat_GM (geometric) : %10.3f\n", Y_hat_GM))
  cat(sprintf("Y_hat_AM (arithmetic): %10.3f\n", Y_hat_AM))
  cat("\n")
  cat(sprintf("Absolute error GM    : %10.3f\n", err_GM))
  cat(sprintf("Absolute error AM    : %10.3f\n", err_AM))
}

Y_hat_GM (geometric) :    476.167
Y_hat_AM (arithmetic):    481.431

Absolute error GM    :      0.000
Absolute error AM    :      5.264

Question 2 – Dataset: airquality

data <- na.omit(airquality[, c("Ozone", "Solar.R", "Wind", "Temp")])

# (a) mat_mult() manual
mat_mult <- function(A,B){
  m <- nrow(A)
  n <- ncol(A)
  p <- ncol(B)
  C <- matrix(0,nrow=m,ncol=p)
  
  for(i in 1:m){
    for(j in 1:p){
      val <- 0
      for(k in 1:n){
        val <- val+A[i,k]*B[k,j]
      }
      C[i,j] <- val
    }
  }
  return(C)
}

# validation
a <- matrix(c(1, 2, 3, 4), nrow = 2)
b <- matrix(c(5, 6, 7, 8), nrow = 2)

c_std    <- a %*% b
c_manual <- mat_mult(a, b)

{
  cat("Standard result (a %*% b):\n")
  print(c_std)
  cat("\nManual result mat_mult(a, b):\n")
  print(c_manual)
  cat("\nElement-wise match:\n")
  print(c_std == c_manual)
}

Standard result (a %*% b):
     [,1] [,2]
[1,]   23   31
[2,]   34   46

Manual result mat_mult(a, b):
     [,1] [,2]
[1,]   23   31
[2,]   34   46

Element-wise match:
     [,1] [,2]
[1,] TRUE TRUE
[2,] TRUE TRUE

# (b) corr_matrix()
corr_matrix <- function(X){
  X <- as.matrix(X)
  n <- nrow(X)
  p <- ncol(X)
  
  # (i) column mean vector X_bar
  ones <- matrix(1,nrow=n,ncol=1)
  X_bar <- (1/n)*mat_mult(t(X),ones)
  
  # (ii) centered matrix X_c and covariance matrix S
  X_c <- X-mat_mult(ones,t(X_bar))
  S <- (1/(n-1))*mat_mult(t(X_c),X_c)
  
  # (iii) D^{-1/2} and correlation matrix R
  D_inv_half <- diag(1/sqrt(diag(S)))
  
  # Return R
  R <- mat_mult(mat_mult(D_inv_half,S),D_inv_half)
  rownames(R) <- colnames(X)
  colnames(R) <- colnames(X)
  return(R)
}

# apply to air quality data and display
R_hat <- corr_matrix(data)
print(round(R_hat, 6))

            Ozone   Solar.R      Wind      Temp
Ozone    1.000000  0.348342 -0.612497  0.698541
Solar.R  0.348342  1.000000 -0.127183  0.294088
Wind    -0.612497 -0.127183  1.000000 -0.497190
Temp     0.698541  0.294088 -0.497190  1.000000

# validation
cat("\nElement-wise match with cor():\n\n")

Element-wise match with cor():

print(round(R_hat, 10) == round(cor(data), 10))

        Ozone Solar.R Wind Temp
Ozone    TRUE    TRUE TRUE TRUE
Solar.R  TRUE    TRUE TRUE TRUE
Wind     TRUE    TRUE TRUE TRUE
Temp     TRUE    TRUE TRUE TRUE

# (c) corr_test_matrix() via trapezoidal integration
corr_test_matrix <- function(X,alpha=0.05){
  X <- as.matrix(X)
  n <- nrow(X)
  p <- ncol(X)
  df <- n-2
  N  <- 10000
  
  # correlation matrix R
  R <- corr_matrix(X)
  
  # initialise p-value matrix
  P <- matrix(0,nrow=p,ncol=p)
  rownames(P) <- colnames(X)
  colnames(P) <- colnames(X)
  
  # p-values for each off-diagonal pair
  for(i in 1:p){
    for(j in 1:p){
      if(i==j){
        
        # set diagonal p-value to 0
        P[i, j] <- 0
        
      } else {
        
        r_ij <- R[i,j]
        
        # (i) test statistic t_ij
        t_ij <- r_ij*sqrt(n-2)/sqrt(1-r_ij^2)
        
        # (ii) trapezoidal integration of dt(t, df) from 0 to |t_ij|
        abs_t <- abs(t_ij)
        h <- abs_t/N
        total <- dt(0,df)+dt(abs_t,df)
        for(k in 1:(N-1)){
          total <- total+2*dt(k*h, df)
        }
        integral <- (h/2)*total
        
        # (ii) p_ij = 1 - 2 * integral
        p_ij <- 1-2*integral
        
        # store result in P[i, j]
        P[i,j] <- p_ij
      }
    }
  }
  return(list(R = R, P = P))
}

# apply to air quality data
results <- corr_test_matrix(data, alpha = 0.05)
{
cat("========== Correlation Matrix (R) ==========\n")
print(round(results$R, 6))

cat("\n========== P-value Matrix (P) ==========\n")
print(round(results$P, 6))
}

========== Correlation Matrix (R) ==========
            Ozone   Solar.R      Wind      Temp
Ozone    1.000000  0.348342 -0.612497  0.698541
Solar.R  0.348342  1.000000 -0.127183  0.294088
Wind    -0.612497 -0.127183  1.000000 -0.497190
Temp     0.698541  0.294088 -0.497190  1.000000

========== P-value Matrix (P) ==========
           Ozone  Solar.R     Wind     Temp
Ozone   0.000000 0.000179 0.000000 0.000000
Solar.R 0.000179 0.000000 0.183452 0.001731
Wind    0.000000 0.183452 0.000000 0.000000
Temp    0.000000 0.001731 0.000000 0.000000

# validation: compare P against cor.test() p-values
cor_mat <- matrix(0, ncol=4, nrow=4)
for (i in 1:4){
  for (j in 1:4){
    if(i==j){
      cor_mat[i,j] <- 0
    } else cor_mat[i,j] <- cor.test(data[,i], data[,j])$p.value
  }
}
rownames(cor_mat) <- colnames(cor_mat) <- names(data)

{
  cat("\nCorrelation Matrix\n")
  print(cor(data))
  cat("\n")
  cat("P-value matrix\n")
  print(cor_mat)
}

Correlation Matrix
             Ozone    Solar.R       Wind       Temp
Ozone    1.0000000  0.3483417 -0.6124966  0.6985414
Solar.R  0.3483417  1.0000000 -0.1271835  0.2940876
Wind    -0.6124966 -0.1271835  1.0000000 -0.4971897
Temp     0.6985414  0.2940876 -0.4971897  1.0000000

P-value matrix
               Ozone      Solar.R         Wind         Temp
Ozone   0.000000e+00 0.0001793109 9.089415e-13 1.552677e-17
Solar.R 1.793109e-04 0.0000000000 1.834520e-01 1.730786e-03
Wind    9.089415e-13 0.1834519762 0.000000e+00 2.841966e-08
Temp    1.552677e-17 0.0017307856 2.841966e-08 0.000000e+00

# verify answer
round(results$P,10) == round(cor_mat,10) 

        Ozone Solar.R  Wind  Temp
Ozone    TRUE    TRUE  TRUE  TRUE
Solar.R  TRUE    TRUE FALSE FALSE
Wind     TRUE   FALSE  TRUE  TRUE
Temp     TRUE   FALSE  TRUE  TRUE