#!/usr/local/bin/R

# Xavier Fernández i Marín
# August 2006
# http://xavier-fim.net

# Some exercises with matrices and matrix algebra


# Simulate independent variables
obs <- 20                    # define a number of observations
X1 <- matrix(rnorm(obs))     # create two random vectors
X2 <- matrix(rnorm(obs))


# Simulate a relationship with the independent variables and Y
alpha <- 10
beta.1 <- 3
beta.2 <- 2
Y <- alpha + beta.1*X1 + beta.2*X2 + rnorm(obs)



# Solve the regression using lm()
model.1 <- lm(Y ~ X1 + X2)
summary(model.1)




# Solve the regression using matrix algebra, step by step
# Least squares estimates of alpha, beta.1 and beta.2 are:
#
# b = [ (X^Transpose * X) ^ Inverse ] * [ X^Transpose * Y ]
#
constant <- matrix(rep(1, obs))        # define the constant as a repetition of 1's
X <- matrix(c(X1, X2, constant), ncol=3)

X.transposed <- t(X)
X.transposed.times.X <- X.transposed %*% X
X.transposed.times.X.inverse <- solve(X.transposed.times.X)

X.transposed.times.Y <- X.transposed %*% Y

B <- X.transposed.times.X.inverse %*% X.transposed.times.Y
B

# Using matrix algebra again, but in one single line
B <- solve(t(X) %*% X )  %*%  (t(X) %*% Y)
B


# Calculate standard errors
df <- obs - (ncol(cbind(X, Y)) - 1)
Y.fitted <- X %*% B
Errors <- Y - Y.fitted
RSS <- t(Errors) %*% Errors
S2 <- RSS/df
S2

V <- S2 kronecker(S2, X.transposed.times.X.inverse)
V <- S2 %x% X.transposed.times.X.inverse # two different ways to calculate
                                         # Kronecker product
V
SE <- sqrt(t(diag(V)))
SE <- t(SE)
SE

# (idem, but in one single line
SE <- sqrt(t(diag(kronecker( (t(Y-(X %*% B)) %*% (Y-(X %*% B))) /
                               (obs - (ncol(cbind(X, Y)) - 1)),
                            solve(t(X) %*% X )))) )

# Calculate t-statistic
t.stat <- B / SE

# Probability of t-statistic
prob.t <- pt(t.stat, df=df, lower.tail=F)*2


# Print the results nicely
results <- data.frame(var=c("alpha", "beta.1", "beta.2"),
                      coef=round(B, 2),
                      std.err=SE,
                      t.val=t.stat,
                      sign=prob.t)
results