Matrix Operations in Mata

* This post demonstrates a few methods for how to input matrices into Mata and how to do some basic matrix operations.  Mata is a matrix programming language so basic matrix operations are extremely easy.

mata
// Let's first build some matrices in Mata
// They can be built directly
A = ( 2 , -1, 5 \ 3, 0 , -1 \ 3, 3 , -1)
A

B1 = (0 , 3 ,-1)
B2 = (3 , 0 , 0)
B3 = (0 , 2 , 0)

// Or through a combination of vectors
B = (B1 \ B2 \ B3)

B

// We can also start with an empty matrix and fill in values

C = J(3,2,4)
// The J command creates a matrix with 3 rows, 2 colums, and with default values of 4

C
// We can replace individual elements once the matrix is created
C[1,2] = 3
C[3,1] = 7

// Or entire submatrices
C[2,] = (2,6)
C

// Now let's see how various matrix operations perform in Mata:

// a. AB
A*B

// b. BA
B*A

// c. A+B
A+B

// d. A'B'
A'*B'

// e. B'A'
B'*A'

// f. A'B
A'*B

// g. (AB)'
(A*B)'

// h. ABC
A*B*C

// i. C'AC
C'*A*C

// j. CAC'
C*A*C'
// This does not exist because C*A cannot be multiplied as C is 3x2 and A is 3x3

// k. trace(C'AC)
trace(C'*A*C)

// l. trace(CAC')
trace(C*A*C')
// Likewise the trace does not exist

end

* Matrices can also be input into Mata from data sets.

clear
set obs 5
gen y = 3
replace y = 6  if _n == 2
replace y = 10 if _n == 3
replace y = 8  if _n == 4
replace y = 2  if _n == 5

mata
// The command st_data retrieves data from stata variables to be used in Mata
y = st_data(. , "y")
y

// Once input, matrices can easily be manipulated
y*y'

// The square of the norm of y
norm(y)^2

// This happens to be identical to:
sum(y:^2)

// The norm is the equclidean norm which is the square square root of the sum of all of the squares of a vector.
// Thus the square of it is just the sum of the squares.

// As should be clear, manipulating matrices in Mata is extremely easy.
// Thus Stata is able to pack a powerful Matrix programming language inside an effective high level user language.

// The largest frustration that I have had with Mata is the relative quality of the documentation.
// I find the documentation of Mata much harder to use than that of Stata (at least in version 11, perhaps version 12 has better documentation).
end

Matrix operations in R

# There are many ways of inputing matrices into R

# Cbind will bind vectors or other matrices together by adding on the columns of one to the other.
A = cbind(c(2,3),c(3,2),c(1,3))

# by the way, the c() function binds a set of elements seperated by commas into a vector.

A

# Rbind does the same as cbind but uses rows instead
B = cbind(c(-1, 1, 4), c(0, 1, 5), c(3, -1, 1))

B

# We can also create a matrix by using the matrix command.
# Data input in this manner is read from a vector into columns.

C = matrix(c(1,4,3,2),nrow = 2, ncol=2)

C

# An array can accompish the same effect as a matrix.

# The largest difference is that an array can take on more than two dimensions.

D = array(c(10,4,5,2), dim=c(2,2))

D

# One need not populate an matrix/array at creation when specifying matrix or array.

E = matrix( NA, nrow = 4, ncol = 4)

E
# We can see the matrix is filled with empty values

# We can replace individual elements by specifying their positions
E[1,1] = 3
E[2,1] = 0
E[3,1] = 4
E[4,1] = -1

# As well as entire columns or rows
E[,2] = c(0,2,3,2)

# Or subsection of the matrix with another matrix
E[,3:4] = cbind(c(2,3,2,1),c(-1,2,1,0))

E

# Sometimes you do not need to specify every element of the matrix if there are common elements
F = matrix(0, nrow=3, ncol=3)

F[1,1] = 4
F[2,2] = 3
F[3,3] = 6

F
# F is a diagnol matrix which is populated first with 0s then filled with the diagnol values.

# Sometimes we start with a data set and want to convert it to a matrix
G = data.frame(score1=c(5,10,-7),score2=c(-1,17,3),score3=c(0,5,2))

G

# Now let's convert it to a matrix
G = data.matrix(G)

G

# Now let's do some matrix operations with the matrices we have defined:

# a. C x D

C
D

C %*% D

# Which is different from element wise multiplication which is

C * D

# b. A x B

A
B

A %*%B

# c. B x A

# d. E x E' = E times the transpose of E

E
E %*% t(E)

# e. F^-1 or F inverse.  We know the inverse of F exists because it is a diagnol matrix with poisitive diagnol non-zero elements.

solve(F)

require(MASS)
ginv(F)

# Neither commands work for this particular application.

# However, it is easy to see that in a diagnol matrix the inverse is just.

# First we will make Finv the same as F
Finv = F

# This is an interesting bit of R functionality

# On the right the diag command is retrieving the diagnol of F and doing a element wise 1/x operation.

# The diag on the left is specifying target values to be replaced.
diag(Finv) = 1/diag(F)

# This might seem a little odd.  Let's try this kind of thing on Z
Z <- Finv

Finv

diag(Z) <- 23
Z
# In this case because 23 is a single number that can be duplicated throughout the diagnol of Z, there is no problem.


# We can check that this is really the inverse of F

Finv %*% F

# Or equally good

F %*% Finv

# Interestingly element by element multiplation yeilds the same result in this example

F * Finv

# f. Rank(D)

D
qr(D)$rank

# F only has rank 1.  This can be seen by dividing col 1 by col 2.

D[,1]/D[,2]

# Both numbers are two meaning item [1,1] = [1,2]*2 and [2,1]=[2,2]*2

# Interestingly this trick works for rows as well:

D[1,]/D[2,]

# g. B + G

B
G

B + G

# h. B - G

B - G

# Addition and subtraction is element by element in matrix notation

# i. Rank(C)

C
qr(C)$rank

# C does not suffer from linearity problems
C[1,]/C[2,]

# j. -G
G
-G

# k. Trace(E)

# Trace is the sum of diagnol elements

E
sum(diag(E))

# l. Rank(F)

# Because F is a diagnol matrix (with no zero diagnol elements) the rank must be equal to the min of the dimensions, 3.
qr(F)$rank

# kF, where k = -7

-7*F

# n. BF

B
F
B %*% F

# o. FB

F %*% B
# With matrices, BF != FB

# p. determinate of C or |C|
C
det(C)

# q. |D|
D
det(D)
# Which looks kind of funny but that is because R usings search algorithms to find the determinate and they are not exact.

# But calculating the determinate of a 2d matrix is easy:
D[1,1]*D[2,2]-D[1,2]*D[2,1]

# r. | CD |
C; D
H = C %*% D
det( H )
# Which also happens to be zero

H[1,1]*H[2,2]-H[1,2]*H[2,1]

# Sorry for the repetitive examples.  I had homework for a class and thought I might as well turn it into a post that someone might find useful.

Write your own System IV estimator in R (and SOLS)

# Specify a variable to hold number of observations
obs = 10000

# Create independent variables
z1 = rnorm(obs)
z2 = rnorm(obs)

# Create error
u1 = rnorm(obs)
u2 = rnorm(obs)

# Create endogenous variable
x1 = .5*z1  -   z2 + u1 + rnorm(obs)
x2 = 1.5*z1 - 2*z2 - u2 + rnorm(obs)
x3 = rnorm(obs)

# Create dependent variable
y1 =  5 + 2*x1 - x2 +   x3 + u1*5 - u2*2
y2 = -5 - 2*x1 + x2 - 2*x3 + u2*5 - u1*1

# First let's attempt a nieve SOLS
# Y = XB
# X'Y = X'XB
# (X'X)^-1(X'Y) = Bhat

X = cbind(x1,x2,x3,1)
Y = cbind(y1,y2)

A = solve(t(X)%*%X)
# solve simply finds the inverse of its argument
B = t(X)%*%Y

SOLS = A%*%B
SOLS

# We can see the coefficient on x is upwards biased.

# Now let's construct the system IV

# Y = XB
# Z'Y = Z'XB
# (Z'X)'Z'Y = (Z'X)'Z'XB
# (X'Z)Z'Y  = (X'Z)Z'XB
# ((X'Z)Z'X)^-1 (X'Z)Z'Y = B
# ((X'Z Z'X)^-1 (X'Z Z'Y) = IVhat
# ((X'Z Z'X)^-1 (X'Z Z'Y) = IVhat
# C * D = IVhat

X = cbind(x1, x2, x3, 1)

# It is important to remember that any explanatory variables that are not instrumented for must be included in the instrument.  Thus x3 in Z.
Z = cbind(x3, z1, z2, 1)

C = solve(t(X)%*%Z%*%t(Z)%*%X)
D = t(X)%*%Z%*%t(Z)

SIV = C%*%D%*%Y
SIV

# System IV seems to be working pretty well
# Alternative non-system IV on the two equations would be

IV1 = C%*%D%*%y1
IV1

IV2 = C%*%D%*%y2
IV2

# We can see that there is no difference in coefficient estimates between seperate IV estimates and SIV.  This, I believe, is because like in the SOLS model if the regressors are the same in both equations and there is no cross equation restrictions then the system produces identical results to the non-system estimates.  This however is not true when estimating standard errors.