#
# Problem is on EPage 203
# Numerical estimates from Example 5.1 on EPage 169 - 170
#
N = 10
m = 2
r = 0.95
r1y = 0.85
r2y = 0.78

S_yy = 1

print('Model X1 and X2: The two predictor case i.e. model is y(x)= beta_1 x_1 + beta_2 x_2 + e')

XtX = matrix( c( 1, r, r, 1 ), nrow=2, ncol=2, byrow=T )
Xty = matrix( c( r1y, r2y ), nrow=2, ncol=1 )
XtX_Inverse = solve(XtX)

beta_hat = solve( XtX, Xty )
print('Model X1 and X2: beta hat=')
print( beta_hat )

# From EPage 126:
#
# REGSS = beta_c^T X_c^T beta_c
#
REGSS = t(beta_hat) %*% XtX %*% beta_hat
REGSS = as.double(REGSS)

# From EPage 126:
#
# RSS = S_yy - REGSS
#
RSS = S_yy - REGSS

# From EPage 46:
#
# s^2 = RSS/(N-p)
#
s2 = RSS/(N-3) # think should be 2 (the mean is assumed zero and we are only estimating beta_1 and beta_2 or two numbers)
print(sprintf('Model X1 and X2: s2= %f', s2))

var_beta_hat = s2 / (1-r^2)
print(sprintf('Model X1 and X2: var_beta_hat= %f', var_beta_hat))

#
# In general var_beta_hat_i = sigma^2 (XtX_Inverse)(i,i)
# See Regression Analysis by Example by Chattergee
#
se_beta_hat = sqrt(var_beta_hat) # the standard error of \hat{beta}
print(sprintf('Model X1 and X2: se_beta_hat= %f', se_beta_hat))


x_sample = c( 0.2, 0.7 ) # what to predict with the models below
#x_sample = c( 0.5, 0.5 )

y_hat = as.double( t(beta_hat) %*% matrix(x_sample, nrow=2, ncol=1)) # our prediction

ts = beta_hat / se_beta_hat
print('Model X1 and X2: ts=')
print(ts)

F = REGSS/(m*s2)
R2 = REGSS / S_yy
print(sprintf('Model X1 and X2: F= %f; R2= %f', F, R2))

# From EPage 130
#
var_prediction = s2 * ( 1 + 1/N + t(x_sample) %*% XtX_Inverse %*% x_sample ) # error in that prediction
se_prediction = sqrt( var_prediction * qf(1-0.05, 1, N-2) )
print(sprintf('Model X1 and X2: y_hat= y(%.1f, %.1f)= %f; se_prediction= %f', x_sample[1], x_sample[2], y_hat, se_prediction))


print('Model X1 only: The single (x_1) predictor case i.e. model is y(x) = beta_1 x_1 + e')

beta_hat = r1y/1

REGSS = beta_hat * 1 * beta_hat

RSS = S_yy - REGSS
s2 = RSS/(N-2)
se_beta_hat = sqrt(s2)
t = beta_hat / se_beta_hat
R2 = REGSS / S_yy
print(sprintf('Model X1 only: beta_1_hat= %f; s2= %f; t= %f; R2= %f', beta_hat, s2, t, R2)) 


x = x_sample[1]
y_hat = beta_hat * x
var_prediction = s2 * ( 1 + 1/N + x * 1 * x ) # error in that prediction
se_prediction = sqrt( var_prediction * qf(1-0.05, 1, N-1) )
print(sprintf('Model X1 only: y_hat= y(%.1f)= %f; se_prediction= %f', x, y_hat, se_prediction))


print('Model X2 only: The single (x_2) predictor case i.e. model is y(x) = beta_2 x_2 + e')

beta_hat = r2y/1

REGSS = beta_hat * 1 * beta_hat
RSS = S_yy - REGSS
s2 = RSS/(N-2)
se_beta_hat = sqrt(s2)
t = beta_hat / se_beta_hat
R2 = REGSS / S_yy

print(sprintf('Model X2 only: beta_2_hat= %f; s2= %f; t= %f; R2= %f', beta_hat, s2, t, R2))

x = x_sample[2]
y_hat = beta_hat * x
var_prediction = s2 * ( 1 + 1/N + x * 1 * x ) # error in that prediction
se_prediction = sqrt( var_prediction * qf(1-0.05, 1, N-1) )

print(sprintf('Model X2 only: y_hat= y(%.1f)= %f; se_prediction= %f', x, y_hat, se_prediction))