> 
> n <- length(Y)
> 
> X0 <- rep(1,n)
> 
> X <- as.matrix(cbind(X0,X))    # Form the X-matrix (n=25 rows, 2 Cols)
> 
> Y <- as.matrix(Y,ncol=1)       # Form the Y-vector (n=25 rows, 1 col)
> 
> # Notes:   t(X) = transpose of X, %*% = matrix multiplication, solve(A) = A^(-1) 
> 
> (XX <- t(X) %*% X)             # Obtain X'X matrix (2 rows, 2 cols)
     X0      X
X0   25   1750
X  1750 142300
> 
> (XY <- t(X) %*% Y)             # Obtain X'Y vector (2 rows, 1 col)
     [,1]
X0   7807
X  617180
> 
> (XXI <- solve(XX))             # Obtain (X'X)^(-1) matrix (2 rows, 2 cols)
             X0             X
X0  0.287474747 -3.535354e-03
X  -0.003535354  5.050505e-05
> 
> (b <- XXI %*% XY)              # Obtain b-vector (2 rows, 1 col)
        [,1]
X0 62.365859
X   3.570202
> 
> Y_hat <- X %*% b               # Obtain the vector of fitted values (n=25 rows, 1 col)
> 
> e <- Y - Y_hat                 # Obtain the vector of residuals (n=25 rows, 1 col)
> 
> print(cbind(Y_hat,e))
          [,1]        [,2]
 [1,] 347.9820  51.0179798
 [2,] 169.4719 -48.4719192
 [3,] 240.8760 -19.8759596
 [4,] 383.6840  -7.6840404
 [5,] 312.2800  48.7200000
 [6,] 276.5780 -52.5779798
 [7,] 490.7901  55.2098990
 [8,] 347.9820   4.0179798
 [9,] 419.3861 -66.3860606
[10,] 240.8760 -83.8759596
[11,] 205.1739 -45.1739394
[12,] 312.2800 -60.2800000
[13,] 383.6840   5.3159596
[14,] 133.7699 -20.7698990
[15,] 455.0881 -20.0880808
[16,] 419.3861   0.6139394
[17,] 169.4719  42.5280808
[18,] 240.8760  27.1240404
[19,] 383.6840  -6.6840404
[20,] 455.0881 -34.0880808
[21,] 169.4719 103.5280808
[22,] 383.6840  84.3159596
[23,] 205.1739  38.8260606
[24,] 347.9820  -5.9820202
[25,] 312.2800  10.7200000
> 
> H <- X %*% XXI %*% t(X)              # Obtain the Hat matrix
> 
> J_n <- matrix(rep(1/n,n^2),ncol=n)   # Obtain the (1/n)J matrix (n=25 rows, n=25 cols)   
> 
> I_n <- diag(n)                       # Obtain the identity matrix (n=25 rows, n=25 cols)  
> 
> (SSTO <- t(Y) %*% (I_n - J_n) %*% Y) # Obtain Total Sum of Squares (SSTO) 
       [,1]
[1,] 307203
> # SSTO can also be computed as:
> # (SSTO <- (t(Y) %*% Y) - (t(Y) %*% (I_n - J_n) %*% Y))
> 
> (SSE <- t(Y) %*% (I_n - H) %*% Y)    # Obtain Error Sum of Squares (SSE)             
         [,1]
[1,] 54825.46
> # SSE can also be computed as:
> # (SSE <- (t(Y) %*% Y) - (t(b) %*% XY))
> 
> (SSR <- t(Y) %*% (H - J_n) %*% Y)    # Obtain Regression Sum of Squares (SSR)             
         [,1]
[1,] 252377.6
> # SSR can also be computed as:
> # (SSR <- (t(b) %*% XY) - (t(Y) %*% J_n %*% Y))
> 
> (MSE <- SSE/(n-2))                   # Obtain MSE = s^2
         [,1]
[1,] 2383.716
> 
> (s2_b <- MSE[1,1] * XXI)             # Obtain s^2{b}, must use MSE[1,1] and * to do scalar multiplication
           X0          X
X0 685.258045 -8.4272774
X   -8.427277  0.1203897
> 
> (X_h <- matrix(c(1,65),ncol=1))      # Create X_h vector, for case where lot size=65
     [,1]
[1,]    1
[2,]   65
> 
> (Y_hat_h <- t(X_h) %*% b)            # Obtain the fitted value when lot size=65
        [,1]
[1,] 294.429
> 
> (s2_yhat_h <- t(X_h) %*% s2_b %*% X_h)        # Obtain s^2{Y_hat_h}
         [,1]
[1,] 98.35837
> 
> (s2_pred <- MSE + (t(X_h) %*% s2_b %*% X_h))  # Obtain s^2{pred}
         [,1]
[1,] 2482.074