Doing Matrix Algebra with Limdep

The purpose of this assignment is to introduce you to some of the capabilities of LIMDEP for doing Matrix Algebra. Below are some exercises which you may enter manually. Note that if you get into trouble you can also cut and paste from the web page into LIMDEP. Table A.1 near the middle of the exercise will need to be entered manually into a lotus worksheet.

/* The reset command clears all of LIMDEP's storage areas */

RESET $

/* Now, let's replicate some materials in Appendix A of Green dealing with
matrix algebra. First, let's experiment with how to create and manipulate
matrices. A matrix is just a rectangular array of numbers with dimension
r x c, where r is the number of rows and c is the number of columns. Let's
create a 3 x 3 symmetrical matrix. Note that in order to see the matrix you
must use the LIST option prior to creating
the matrix. */

MATRIX ; LIST ; A=[1,3,7/3,5,2/7,2,1] $

/* Identity matrices are created using a special command within LIMDEP. Let's
create a 3 x 3 identity matrix. */

MATRIX ; LIST ; I3=IDEN(3) $

/* Manipulation of matrices is easy within LIMDEP. One important operation
is transposition, which means to switch the columns with the rows of the matrix.
Let's create a new matrix and transpose it. */

MATRIX ; LIST ; A=[1,2,3/5,1,5/6,4,5/3,1,4] $
MATRIX ; LIST ; AT=A' $

/* Matrix addition and subtraction means to add or subtract all corresponding
elements of two matrices. Addition and subtraction require matrices of the
same dimensions. Let's create a new matrix C, and add it to A */

MATRIX ; LIST ; C=[1,1,1/2,2,2/3,3,3/4,4,4] $
MATRIX ; LIST ; AplusC=A+C $

/* Matrix multiplication involves taking the inner or dot products of rows and
columns. Let's begin by taking the inner product of two vectors to obtain a scalar. */

MATRIX ; LIST ; A=[1,3,4] ; C=[3/8/2] $
MATRIX ; LIST ; AC=A*C $

/* Matrices of higher dimension require conformability. Conformability means
that the number of columns of the first matrix equals the number of rows for
the second matrix. */

MATRIX ; LIST ; A=[1,3,2/4,5,-1] ; C=[2,4/1,6/0,5] $
MATRIX ; LIST ; AC=A*C $

/* Matrix multiplication is generally not commutative, but it can be for certain
matrices. */

MATRIX ; LIST ; CA=C*A $

/* Matrix multiplication is just a linear combination of the columns of the first
matrix with the rows of the second matrix. For example, take the matrix AC above.
We can break C down into column 1 and column 2, called C1 and C2 and perform
the multiplications separately. */

MATRIX ; LIST ; C1=[2/1/0] ; C2=[4/6/5] $

/*Then performing the multiplications we have two column vectors which comprise AC. */

MATRIX ; LIST ; AC1=A*C1 ; AC2=A*C2 $
MATRIX ; LIST ; AC=A*C $

/* The identity matrix performs the function of the number 1 in matrix algebra.
Note that the identity matrix is an example of an idempotent matrix, where
the product is always commutative. An idempotent matrix is equal to its square.*/

MATRIX ; LIST ; A=[4,2,1/2,6,1/1,1,0] ; I=IDEN(3)$
MATRIX ; LIST ; AI=A*I ; IA=I*A $

/* Matrix products are associative and distributive */

MATRIX ; LIST ; A=[1,2/3,4] ; C=[2,3/4,5] ; D=[3,4/5,6] $
MATRIX ; LIST ; AC=A*C ; CD=C*D ; ACtimesD=AC*D ; AtimesCD=A*CD $
MATRIX ; LIST ; CplusD=C+D ; AC=A*C ; AD=A*D ; ACplusD=A*CplusD ; ACplusAD=AC+AD $

/*A useful bit of information is that (AC)'=C'A', and (ACD)'=D'C'A'. This
generalizes to any number of products. */

MATRIX ; LIST; ACtransp=AC' ; CtrAtr=C'*A' $

/* Some other useful bits of information are how to compute n, sums of columns,
means, sums of squares, and cross products */

MATRIX ; LIST ; ONES=INIT(5,1,1) $

/*Compute N */

MATRIX ; LIST ; NUM=ONES'ONES $

/* Compute sum of the elements of a vector x */

MATRIX ; LIST ; x=[1/2/3/4/5] $
MATRIX ; LIST ; SumX=ONES'x $

/* Compute sum of all ax */

MATRIX ; LIST ; a=.5 ; sumax=a*ones'x $

/* Compute the sum of squares for x */

MATRIX ; LIST ; XSS=x'x $

/* Compute the cross product of x and another vector y */

MATRIX ; LIST ; y=[2/3/4/5/6] $
MATRIX ; LIST ; xdoty=x'y $

/* Now, let's input the data matrix from Table A.1 of Greene, and create
a matrix using the namelist command. */

RESET $
READ ; FORMAT=WKS ; NAMES ; FILE=TBLA1.WK1 $
DSTAT ; RHS=* $
NAMELIST ; X=YEAR, CONSUME, GNP $

/* Now let's take the inner product of the matrix X */

MATRIX ; LIST ; XTX=X'X $

/*Note that on the diagonal we have the sums of squares
for year, consume, and gnp; off the diagonal we have the cross
products of year, consume, and gnp. For example, row 2
and column 3 marks the cross product of consume and gnp,
as do row 3 and column 2. */

/* A useful idempotent matrix, M0, is the one that is used
to transform variables to deviation units from their
mean. Let's create that matrix for the data in Table 2.1 . */

MATRIX ; LIST ; ONES=INIT(9,1,1) $
MATRIX ; LIST ; NUM=ONES'ONES $
MATRIX ; LIST ; ONEOVERN={1/NUM}*ONES*ONES' $
MATRIX ; LIST ; I=IDEN(9) $
MATRIX ; LIST ; M0=I-ONEOVERN $

/* Now lets do the transformation of X to deviations */

MATRIX ; LIST ; XDEV=M0*X $

/* M0 is idempotent, which means that it is also equal
to its square. It is also equal to M0'M0. */

MATRIX ; LIST ; M0M0=M0*M0 $
MATRIX ; LIST ; M0TM0=M0'M0 $

/* M0 is also useful in simplifying operations. In particular,
suppose we want X'X in deviation form, but don't want to
convert the actual data, then do the following. */

MATRIX ; LIST ; XTXDEV=X'*M0*X $

/* Rank of a matrix. The rank of a matrix is the number of
linearly independent columns or rows. The rank of a matrix
cannot exceed the smallest of the number of columns or rows.
It may, however, be less if two or more columns or rows are
linearly dependent. */

MATRIX ; LIST ; C=[1,2,3/5,1,5/6,4,5/3,1,4] $
MATRIX ; LIST ; CRANK=RANK(C) $
MATRIX ; LIST ; CTRANK=RANK(C') $

/* Determinant of a square matrix. The determinant of a square
matrix is a number that reflects the degree of linear dependence
among the vectors (columns) of the matrix. A determinant of zero
says that one or more rows are perfectly linearly dependent.
A determinant is simply the area or volume of the n
dimensional space defined by the vectors of a matrix. Multiplying
a matrix by a scalar is equivalent to multiplying the volume
of the box by the scalar, which affects the determinant in
scalar fashion. */

MATRIX ; LIST ; A=[4,2/1,3] $
MATRIX ; LIST ; DETA=DTRM(A) $

/* Inverse of a matrix. The inverse of a matrix is a matrix
such that Ainverse*A=I. Procedures for calculating matrix
inverses are somewhat technical, but easily implemented by hand
for up to about a 5x5 matrix. LIMDEP calculates matrix inverses
quite readily. */

MATRIX ; LIST ; AINV=<A> $

/* Check that it is an inverse. Note that due to rounding
error and imprecision, the result may not be perfectly an
identity matrix. */

MATRIX ; LIST ; AAINV=A*AINV $

/* Properties of matrix inverses. The determinant of the
inverse equals 1 over the determinant of the matrix */

MATRIX ; LIST ; AINVDET=DTRM(AINV) ; ADET=DTRM(A) $

/* When all inverses exist ACDinverse=Dinv*Cinv*Ainv */

MATRIX ; LIST ; C=[1,3/2,4] ; D=[5,6/2,1] $
MATRIX ; LIST ; ACD=A*C*D $
MATRIX ; LIST ; ACDINV=<ACD> ; ACDINV2=<D>*<C>*<A> $

/* You can also create partitioned matrices within LIMDEP */

MATRIX ; LIST ; C1=[1,2,3,4/3,5,7,9] ; C2=[0,0,0,0/0,0,0,0]
; C3=[0,0,0,0/0,0,0,0] ; C4=[2,4,6,8/3,5,7,9] $
MATRIX ; LIST ; C1C2C3C4=[C1,C2/C3,C4] $

/* Using the Kronecker product is an easy way to create
certain partitioned matrices such as those required for
panel data. *

MATRIX ; LIST ; I=IDEN(2) ; Z=[1,2/3,4] $
MATRIX ; LIST ; ZI=KRON(I,Z) $

/* Characteristic vectors and roots of a matrix, also called
eigenvectors and eigenvalues, may be easily computed */

MATRIX ; LIST ; A=[5,1/1,4] $
MATRIX ; LIST ; AEIGVEC=CVEC(A) $
MATRIX ; LIST ; I=AEIGVEC*AEIGVEC' $
MATRIX ; LIST ; AEIGVAL=CXRT(A) $

/* We can diagonalize a matrix by the following transformation
where we are using the eigenvectors. */

MATRIX ; LIST ; LAMBDA=AEIGVEC'*A*AEIGVEC $

/* We can retrieve the original matrix through the following,
which is sometimes called the spectral decomposition. This is just
a factorization of the matrix A */

MATRIX ; LIST ; A=AEIGVEC*LAMBDA*AEIGVEC' $

/* Matrix diagonalization is useful in dealing with the autocorrelation
problem in regression. It can also be useful where you need powers of
a matrix, as in Markov Matrices. */

CALC ; LIST ; POWER=2 $
MATRIX ; LIST ; LAMBDAP=LAMBDA^POWER $
MATRIX ; LIST ; AA=AEIGVEC*LAMBDAP*AEIGVEC' ; AA2=A^2 $

/*Note that the determinant of a matrix is the product of the
eigenvalues. If any of the eigenvalues are 0, then the matrix is
singular. Also, the determinant of a symmetrical matrix equals
the determinant of it diagonalized form. */

MATRIX ; LIST ; EIGVAL1=PART(AEIGVAL,1,1,1,1) ; EIGVAL2=PART(AEIGVAL,2,2,1,1) $
MATRIX ; LIST ; DETA=DTRM(A) ; DETA2=EIGVAL1*EIGVAL2 ; DETA3=DTRM(LAMBDA) $

/* The rank of the product of multiple matrices is always less than or
equal to the rank of the matrix with the lowest rank. */

MATRIX ; LIST ; A=[1,2/8,7] ; C=[0,1/5,7] ; D=[2,0,4,4,1/3,1,1,2,2] $
MATRIX ; LIST ; ACD=A*C*D $

MATRIX ; LIST ; RANKACD=RANK(ACD) $

/*Also, for a symmetric matrix, the number of non-zero eigenvalues is the
rank of the matrix. This is useful since rank(A)=rank(A'A), because rank(A'A)
is always symmetrical. This means that you can obtain the rank of any matrix
by taking A'A and finding the number of non-zero eigenvalues. */

MATRIX ; LIST ; XACD=ACD'ACD ; EIGXACD=CXRT(XACD) ; RANK(XACD) $

/* The eigenvalues for the inverse of a matrix are equal to 1 over the eigenvalues
for the matrix inverted */

MATRIX ; LIST ; A=[1,2/8,7] ; AINV=<A> $
MATRIX ; LIST ; AEIG=CXRT(A) ; AINVEIG=CXRT(AINV) $

/* Condition Number of a Matrix: the condition number of a matrix is a measure
of the degree of linear dependence between rows/columns. Let's compute the condition
number for a matrix containing YEAR,CONSUMPTION,GNP,and INTEREST from Table A.1. First,
we norm the data, and then we calculate the condition number.*/

NAMELIST; X=YEAR,CONSUME,GNP,INTEREST $
MATRIX ; LIST ; XTX=X'X $
MATRIX ; LIST ; D=DIAG(XTX) $
MATRIX ; LIST ; D=SQRT(D) $
MATRIX ; LIST ; V=<D>*X'*X*<D> $
MATRIX ; LIST ; EIGVALS=CXRT(V) ; MAXEIG=PART(EIGVALS,1,1,1,1)
; MINEIG=PART(EIGVALS,4,4,1,1) $
MATRIX ; LIST ; CONDNUMR=MAXEIG*DIRI(MINEIG) ; CONDNUMR=SQRT(CONDNUMR) $

/* Trace of a matrix: the trace of a square matrix is the sum of the diagonal elements.
It is used primarily in proofs, but is very easy to compute. */

MATRIX ; LIST ; A=[1,2,3/4,5,6/7,8,9] $
MATRIX ; LIST ; ATRACE=TRCE(A) $

/* Factoring a Matrix: A symmetrical positive definite matrix can be factored
in several ways. We have already noted one factorization of a symmetrical matrix,
namely the 1) Spectral Decomposition, A=C*Lambda*C' where Lambda is a diagonal
matrix containing the eigenvalues on the diagonal and C is a matrix containing the
eigenvectors. Other useful factorizations are 2) Choleski's decomposition, where
A can be written as the product of lower triangular and upper triangular matrices,
3) the Singular Value Decomposition, and 4) the QR factorization. LIMDEP WILL
AUTOMATICALLY COMPUTE CHOLESKI'S FACTORIZATION, but the others have to be done manually.
Procedures for factoring a matrix can be found in any elementary text on linear algebra.*/

MATRIX ; LIST ; A=[2,1,3/1,3,4/3,4,9] $
MATRIX ; LIST ; L=CHOL(A) $
MATRIX ; LIST ; LLT=L*L' $