1 #!/usr/local/bin/R
 2 
 3 # Xavier Fernández i Marín
 4 # August 2006
 5 # http://xavier-fim.net
 6 
 7 # Some exercises with matrices and matrix algebra
 8 
 9 
10 # Simulate independent variables
11 obs <- 20                    # define a number of observations
12 X1 <- matrix(rnorm(obs))     # create two random vectors
13 X2 <- matrix(rnorm(obs))
14 
15 
16 # Simulate a relationship with the independent variables and Y
17 alpha <- 10
18 beta.1 <- 3
19 beta.2 <- 2
20 Y <- alpha + beta.1*X1 + beta.2*X2 + rnorm(obs)
21 
22 
23 
24 # Solve the regression using lm()
25 model.1 <- lm(Y ~ X1 + X2)
26 summary(model.1)
27 
28 
29 
30 
31 # Solve the regression using matrix algebra, step by step
32 # Least squares estimates of alpha, beta.1 and beta.2 are:
33 #
34 # b = [ (X^Transpose * X) ^ Inverse ] * [ X^Transpose * Y ]
35 #
36 constant <- matrix(rep(1, obs))        # define the constant as a repetition of 1's
37 X <- matrix(c(X1, X2, constant), ncol=3)
38 
39 X.transposed <- t(X)
40 X.transposed.times.X <- X.transposed %*% X
41 X.transposed.times.X.inverse <- solve(X.transposed.times.X)
42 
43 X.transposed.times.Y <- X.transposed %*% Y
44 
45 B <- X.transposed.times.X.inverse %*% X.transposed.times.Y
46 B
47 
48 # Using matrix algebra again, but in one single line
49 B <- solve(t(X) %*% X )  %*%  (t(X) %*% Y)
50 B
51 
52 
53 # Calculate standard errors
54 df <- obs - (ncol(cbind(X, Y)) - 1)
55 Y.fitted <- X %*% B
56 Errors <- Y - Y.fitted
57 RSS <- t(Errors) %*% Errors
58 S2 <- RSS/df
59 S2
60 
61 V <- S2 kronecker(S2, X.transposed.times.X.inverse)
62 V <- S2 %x% X.transposed.times.X.inverse # two different ways to calculate
63                                          # Kronecker product
64 V
65 SE <- sqrt(t(diag(V)))
66 SE <- t(SE)
67 SE
68 
69 # (idem, but in one single line
70 SE <- sqrt(t(diag(kronecker( (t(Y-(X %*% B)) %*% (Y-(X %*% B))) /
71                                (obs - (ncol(cbind(X, Y)) - 1)),
72                             solve(t(X) %*% X )))) )
73 
74 # Calculate t-statistic
75 t.stat <- B / SE
76 
77 # Probability of t-statistic
78 prob.t <- pt(t.stat, df=df, lower.tail=F)*2
79 
80 
81 # Print the results nicely
82 results <- data.frame(var=c("alpha", "beta.1", "beta.2"),
83                       coef=round(B, 2),
84                       std.err=SE,
85                       t.val=t.stat,
86                       sign=prob.t)
87 results