You can use Mata interactively.
A = 1,2\3,4 // Input a 2 × 2 matrix
b = -1\3 // Input a 2 × 1 column vector
A'b // Compute A'b
A[2,1] // Access the (2,1) element of A
A[,2] // Access the second column of A
B = diag(b) // Create a 2 × 2 matrix with the elements of b on the diagonal
C = A,B // Combine matrices A and B by column
D = A\B // Or by row
I = I(3) // Create a 3 × 3 identity matrix I
trace(A) // Compute the trace of A
det(A) // And the determinant of A
G*H // Multiply G (n × k) by H (k × k)
G:*g // Multiply each column of G (n × k) by g (n × 1)
st_data(.,.) // Load Stata dataset into Mata
st_data((1,5\7,9), "mpg") // Observations 1 to 5 and 7 to 9 of variable mpg
Or you can write a simple program for regression.
mata
st_view(y=., ., "mpg")
st_view(X=., ., ("weight", "foreign", "cons"))
b = invsym(X'X)*X'y
b
e = y - X*b
n = rows(X)
k = cols(X)
s2= (e'e)/(n-k)
V = s2*invsym(X'X)
V
se = sqrt(diagonal(V))
(b, se)
Or build a regression system using structures or classes.
class linreg
{
public:
void setup()
real scalar N()
real scalar k()
real scalar K()
real colvector b_X()
real scalar b_c()
real colvector b()
real matrix V()
real scalar s2()
real scalar ee()
real colvector yhat()
protected:
// --------------------------------------- inputs
pointer(real colvector) scalar y
pointer(real matrix) scalar X
real scalar cons
// --------------------------------------- derived
real matrix XX // X'X; K x K
real colvector Xy // X'y; K x K
real matrix XXinv // (X'X)^(-1); K x K
real colvector b // coefficients; K x 1
real matrix V // variance matrix (VCE); K x K
real scalar s2 // variance of e
real scalar yy // y'y
real scalar ee // e'e, error sum of squares
real scalar K_adj
real scalar k_adj
/*
Deviation from mean notation:
S = X - mean(X) X as deviation from mean
t = y - mean(y) y as deviation from mean
*/
real matrix SS // S'S; k x k
real colvector St // S't; K x 1 (sic)
real matrix SSinv // (S'S)^(-1); K x K (sic)
real rowvector mean_X // means of rows of X; 1 x k
real scalar mean_y // mean of row of y; 1 x 1
void init()
real matrix XX()
real colvector Xy()
real scalar yy()
real matrix XXinv()
real matrix SSinv()
void calc_SSinv_cons()
real matrix SS()
real colvector St()
real rowvector mean_X()
real scalar mean_y()
real scalar K_adj()
real scalar k_adj()
}
void linreg::init()
{
XX = .z
Xy = .z
XXinv = .z
b = .z
V = .z
s2 = .z
yy = .z
ee = .z
SS = .z
St = .z
mean_X = .z
mean_y = .z
SSinv = .z
K_adj = .z
k_adj = .z
}
/* -------------------------------------------------- entry points --- */
void linreg::setup( real colvector user_y,
real matrix user_X,
real scalar user_cons )
{
assert(rows(y)==rows(X))
y = &user_y
X = &user_X
cons = (user_cons!=0)
init() // initialize derived vars
}
real scalar linreg::N() return(rows(*X))
real scalar linreg::k() return(cols(*X))
real scalar linreg::K() return(cols(*X) + cons)
/*
Throughout this code, distinguish between k and K:
k = # of independent variables EXCLUDING intercept
K = # of independent variables INCLUDING intercept
In models without an intercept, K = k
In models with an intercept, K = k + 1
Note: *X is k x k, X does not include vector of 1s, and
b is K x 1, coefficient include intercept if there is one.
To code the equivalent of X*b, you extract the first k elements
of b and add the intercept separately. To make this easier
think of
b = (b_X, b_c)
b_X() returns b_X: 1 x k
b_c() returns b_c: 1 x 1 containing b[K] or 0.
Thus, X*b can be calculated as (*X)*b_X() :+ b_c()
in all cases.
*/
real colvector linreg::b_X()
{
if (cons) {
if (k()) return( b()[| 1\k() |])
else return( J(0, 1, .))
}
else return(b())
}
real scalar linreg::b_c()
{
return(cons ? b()[K()] : 0)
}
real colvector linreg::b()
{
if (b == .z) {
if (cons) {
b = SSinv() * St()
b[K()] = mean_y() - mean_X()*b_X()
}
else b = XXinv() * Xy()
}
return(b)
}
real matrix linreg::V()
{
if (V == .z) {
V = s2() * (cons ? SSinv() : XXinv())
}
return(V)
}
real scalar linreg::s2()
{
if (s2 == .z) {
s2 = ee() / ( N() - K_adj() )
}
return(s2)
}
real scalar linreg::ee()
{
real colvector e
if (ee == .z) {
e = *y - yhat()
ee = quadcross(e, e)
}
return(ee)
}
real colvector linreg::yhat()
{
return( (*X)*b_X() :+ b_c() )
}
/* ------------------------------------------internal subroutines --- */
real matrix linreg::XX()
{
if (XX == .z) {
XX = quadcross(*X, cons, *X, cons)
}
return(XX)
}
real colvector linreg::Xy()
{
if (Xy == .z) {
Xy = quadcross(*X, cons, *y, 0)
}
return(Xy)
}
real scalar linreg::yy()
{
if (yy == .z) {
yy = quadcross(*y, *y)
}
return(yy)
}
real matrix linreg::XXinv()
{
if (XXinv == .z) {
XXinv = invsym(XX())
}
return(XXinv)
}
real matrix linreg::SSinv()
{
if (SSinv == .z) {
if (cons) calc_SSinv_cons()
}
return(SSinv)
}
void linreg::calc_SSinv_cons()
{
real scalar K, k, i
real rowvector v
real matrix kxk, rowK, colK
K = K() // # X vars including constant
k = k() // # x vars excluding constant
/*
Calculate inverse of S'S = SS().
This routine called only if cons, meaning K = k + 1.
SS() is k x k -- it excludes constant.
We want the full K x K inverse.
*/
// ---------------------------------------- SSinv 1 x 1 case ---
if (K==1) {
SSinv = 1/N()
return
}
// ----------------------------------- SSinv 2 x 2 or larger ---
SSinv = J(K, K, .) // matrix now K x K
// Range subscripts:
kxk = (1,1 \ k,k) // k x k submatrix rows & cols 1..k
rowK = (K,1 \ K,k) // 1 x k submatrix; row K
colK = (1,K \ k,K) // k x 1 submatrix; col K
SSinv[|kxk |] = invsym(SS())
SSinv[|rowK|] = v = -(mean_X() * SSinv[|kxk|])
SSinv[|colK|] = v'
// SSinv[ K,K ] = 1/N() - mean_X()*SSinv[|colK|]
SSinv[ K,K ] = 1/N() - mean_X()*v'
}
real matrix linreg::SS()
{
if (SS == .z) {
if (cons) {
SS = quadcrossdev(*X, 0, mean_X(),
*X, 0, mean_X() )
}
}
return(SS)
}
real colvector linreg::St()
{
if (St == .z) {
St = quadcrossdev(*X, 0, mean_X(),
*y, 0, mean_y()) \ 0
}
return(St)
}
real rowvector linreg::mean_X()
{
if (mean_X == .z) mean_X = mean(*X)
return(mean_X)
}
real scalar linreg::mean_y()
{
if (mean_y == .z) mean_y = mean(*y)
return(mean_y)
}
real scalar linreg::K_adj()
{
if (K_adj == .z) {
if (cons) K_adj = colsum(diagonal(SSinv()) :!= 0)
else K_adj = colsum(diagonal(XXinv()) :!= 0)
}
return(K_adj)
}
real scalar linreg::k_adj()
{
if (k_adj == .z) {
k_adj = k() - (K()-K_adj())
}
return(k_adj)
}
end
Resources
There are many great resources for learning Mata.
The official documentation:
We recommend starting with the intro sections:
The Mata Book: A Book for Serious
Programmers and Those Who Want to Be by William Gould
The Mata Matters articles published in the Stata Journal:
Programming an estimation command in Stata, a series of posts on the Stata Blog.
An Introduction to Mata
from the Social Science Computing Cooperative at the University of Madison–Wisconsin
Learn more about all of Mata's features.
