\chapter{Matrix manipulation}
\label{chap:matrices}
Together with the other two basic types of data (series and scalars),
\app{gretl} offers a quite comprehensive array of matrix methods. This
chapter illustrates the peculiarities of matrix syntax and discusses
briefly some of the more complex matrix functions. For a full listing
of matrix functions and a comprehensive account of their syntax,
please refer to the \GCR.
\section{Creating matrices}
\label{sec:matrix-create}
Matrices can be created using any of these methods:
\begin{enumerate}
\item By direct specification of the scalar values that compose the
matrix --- in numerical form, by reference to pre-existing
scalar variables, or using computed values.
\item By providing a list of data series.
\item By providing a \textit{named list} of series.
\item Using a formula of the same general type that is used
with the \texttt{genr} command, whereby a new matrix is defined
in terms of existing matrices and/or scalars, or via some
special functions.
\end{enumerate}
To specify a matrix \textit{directly in terms of scalars}, the syntax
is, for example:
\begin{code}
matrix A = { 1, 2, 3 ; 4, 5, 6 }
\end{code}
The matrix is defined by rows; the elements on each row are separated
by commas and the rows are separated by semi-colons. The whole
expression must be wrapped in braces. Spaces within the braces are
not significant. The above expression defines a $2\times3$ matrix.
Each element should be a numerical value, the name of a scalar
variable, or an expression that evaluates to a scalar. Directly after
the closing brace you can append a single quote (\texttt{'}) to obtain
the transpose.
To specify a matrix \textit{in terms of data series} the syntax is,
for example,
%
\begin{code}
matrix A = { x1, x2, x3 }
\end{code}
%
where the names of the variables are separated by commas. Besides
names of existing variables, you can use expressions that evaluate to
a series. For example, given a series \texttt{x} you could do
%
\begin{code}
matrix A = { x, x^2 }
\end{code}
%
Each variable occupies a column (and there can only be one variable
per column). You cannot use the semi-colon as a row separator in this
case: if you want the series arranged in rows, append the transpose
symbol. The range of data values included in the matrix depends on
the current setting of the sample range.
\tip{While gretl's built-in statistical functions for data series are
capable of handling missing values, the matrix arithmetic functions
are not. When you build a matrix from series that include missing
values, observations for which at least one series has a missing
value are skipped.}
Instead of giving an explicit list of variables, you may instead
provide the \textit{name of a saved list} (see
Chapter~\ref{chap-persist}), as in
%
\begin{code}
list xlist = x1 x2 x3
matrix A = { xlist }
\end{code}
%
When you provide a named list, the data series are by default placed
in columns, as is natural in an econometric context: if you want them
in rows, append the transpose symbol.
As a special case of constructing a matrix from a list of variables,
you can say
%
\begin{code}
matrix A = { dataset }
\end{code}
%
This builds a matrix using all the series in the current dataset,
apart from the constant (variable 0). When this dummy list is used, it
must be the sole element in the matrix definition \texttt{\{...\}}. You
can, however, create a matrix that includes the constant along with
all other variables using horizontal concatenation (see below), as in
%
\begin{code}
matrix A = {const}~{dataset}
\end{code}
%
The syntax
%
\begin{code}
matrix A = {}
\end{code}
%
creates an empty matrix --- a matrix with zero rows and zero columns.
See section \ref{sec:emptymatrix} for a discussion of this
object.
\tip{Names of matrices must satisfy the same requirements as names of
gretl variables in general: the name can be no longer than 15
characters, must start with a letter, and must be composed of
nothing but letters, numbers and the underscore character.}
\section{Empty matrices}
\label{sec:emptymatrix}
The main purpose of the concept of an empty matrix is to enable the
user to define a starting point for subsequent concatenation
operations. For instance, if \texttt{X} is an already defined
matrix of any size, the commands
%
\begin{code}
matrix A = {}
matrix B = A ~ X
\end{code}
%
result in a matrix \texttt{B} identical to \texttt{X}.
From an algebraic point of view, one can make sense of the idea of an
empty matrix in terms of vector spaces: if a matrix is an ordered set
of vectors, then \verb|A={}| is the empty set. As a consequence,
operations involving addition and multiplications don't have any clear
meaning (arguably, they have none at all), but operations involving
the cardinality of this set (that is, the dimension of the space
spanned by \texttt{A}) are meaningful.
Legal operations on empty matrices are listed in Table
\ref{tab:empty-matrix-funcs}. (All other matrix operations generate
an error when an empty matrix is given as an argument.) In line with
the above interpretation, some matrix functions return an empty matrix
under certain conditions: the functions \texttt{diag, vec, vech,
unvech} when the arguments is an empty matrix; the functions
\texttt{I, ones, zeros, mnormal, muniform} when one or more of the
arguments is 0; and the function \texttt{nullspace} when its argument
has full column rank.
\begin{table}[htbp]
\centering
\begin{tabular}{lc}
\textit{Function} & \textit{Return value} \\ [4pt]
\texttt{A', transp(A)} & \texttt{A} \\
\texttt{rows(A)} & 0 \\
\texttt{cols(A)} & 0 \\
\texttt{rank(A)} & 0 \\
\texttt{det(A)} & \texttt{NA} \\
\texttt{ldet(A)} & \texttt{NA} \\
\texttt{tr(A)} & \texttt{NA} \\
\texttt{onenorm(A)} & \texttt{NA} \\
\texttt{infnorm(A)} & \texttt{NA} \\
\texttt{rcond(A)} & \texttt{NA} \\
\end{tabular}
\caption{Valid functions on an empty matrix, \texttt{A}}
\label{tab:empty-matrix-funcs}
\end{table}
\section{Selecting sub-matrices}
\label{sec:matrix-sub}
You can select sub-matrices of a given matrix using the syntax
\hspace{1em} \texttt{A[}\textsl{rows},\textsl{cols}\texttt{]}
where \textsl{rows} can take any of these forms:
\begin{center}
\begin{tabular}{lll}
1. & empty & selects all rows \\
2. & a single integer & selects the single specified row \\
3. & two integers separated by a colon & selects a range of rows \\
4. & the name of a matrix & selects the specified rows \\
\end{tabular}
\end{center}
With regard to option 2, the integer value can be given numerically,
as the name of an existing scalar variable, or as an expression that
evaluates to a scalar. With the option 4, the index matrix given
in the \textsl{rows} field must be either $p\times 1$ or $1\times p$,
and should contain integer values in the range 1 to $n$, where $n$ is
the number of rows in the matrix from which the selection is to be
made.
The \textsl{cols} specification works in the same way, \textit{mutatis
mutandis}. Here are some examples.
%
\begin{code}
matrix B = A[1,]
matrix B = A[2:3,3:5]
matrix B = A[2,2]
matrix idx = { 1, 2, 6 }
matrix B = A[idx,]
\end{code}
%
The first example selects row 1 from matrix \texttt{A}; the second
selects a $2\times 3$ submatrix; the third selects a scalar; and
the fourth selects rows 1, 2, and 6 from matrix \texttt{A}.
In addition there is a pre-defined index specification, \texttt{diag},
which selects the principal diagonal of a square matrix, as in
\texttt{B[diag]}, where \texttt{B} is square.
You can use selections of this sort on either the right-hand side of
a matrix-generating formula or the left. Here is an example of use of
a selection on the right, to extract a $2\times 2$ submatrix $B$ from a
$3\times 3$ matrix $A$:
%
\begin{code}
matrix A = { 1, 2, 3; 4, 5, 6; 7, 8, 9 }
matrix B = A[1:2,2:3]
\end{code}
%
And here are examples of selection on the left. The second line below
writes a $2\times 2$ identity matrix into the bottom right corner of the
$3\times 3$ matrix $A$. The fourth line replaces the diagonal of $A$
with 1s.
%
\begin{code}
matrix A = { 1, 2, 3; 4, 5, 6; 7, 8, 9 }
matrix A[2:3,2:3] = I(2)
matrix d = { 1, 1, 1 }
matrix A[diag] = d
\end{code}
\section{Matrix operators}
\label{matrix-op}
The following binary operators are available for matrices:
\begin{center}
\begin{tabular}{ll}
\texttt{+} & addition \\
\texttt{-} & subtraction \\
\texttt{*} & ordinary matrix multiplication \\
\texttt{'} & pre-multiplication by transpose \\
\texttt{/} & matrix ``division'' (see below) \\
\verb+~+ & column-wise concatenation \\
\verb+|+ & row-wise concatenation \\
\texttt{**} & Kronecker product \\
\texttt{=} & test for equality
\end{tabular}
\end{center}
In addition, the following operators (``dot'' operators) apply on an
element-by-element basis:
\begin{center}
\begin{tabular}{cccccccc}
\texttt{.+} & \texttt{.-} &
\texttt{.*} & \texttt{./} & \verb+.^+ &
\texttt{.=} & \texttt{.>} & \texttt{.<}
\end{tabular}
\end{center}
Here are explanations of the less obvious cases.
For matrix addition and subtraction, in general the two matrices have
to be of the same dimensions but an exception to this rule is granted
if one of the operands is a $1\times 1$ matrix or scalar. The scalar
is implicitly promoted to the status of a matrix of the correct
dimensions, all of whose elements are equal to the given scalar value.
For example, if $A$ is an $m \times n$ matrix and $k$ a scalar, then
the commands
%
\begin{code}
matrix C = A + k
matrix D = A - k
\end{code}
%
both produce $m \times n$ matrices, with elements $c_{ij} =
a_{ij} + k$ and $d_{ij} = a_{ij} - k$ respectively.
By ``pre-multiplication by transpose'' we mean, for example, that
%
\begin{code}
matrix C = X'Y
\end{code}
%
produces the product of $X$-transpose and $Y$. In effect,
the expression \texttt{X'Y} is shorthand for \texttt{X'*Y}
(which is also valid).
In matrix ``division'', the statement
%
\begin{code}
matrix C = A/B
\end{code}
%
is interpreted as a request to find the matrix $C$ that solves $BC=A$.
If $B$ is a square matrix, this is treated as equivalent to $B^{-1}A$,
which fails if $B$ is singular; the numerical method employed here is
the LU decomposition. If $B$ is a $T \times k$ matrix with
$T > k$, then $C$ is the least-squares solution, $C = (B'B)^{-1}B'A$,
which fails if $B'B$ is singular; the numerical method employed here is
the QR decomposition. Otherwise, the operation necessarily fails.
In ``dot'' operations a binary operation is applied element by
element; the result of this operation is obvious if the matrices are
of the same size. However, there are several other cases where such
operators may be applied. For example, if we write
%
\begin{code}
matrix C = A .- B
\end{code}
%
then the result $C$ depends on the dimensions of $A$ and $B$. Let $A$
be an $m \times n$ matrix and let $B$ be $p \times q$; the result is
as follows:
\begin{center}
\begin{tabular}{ll}
\textit{Case} & \textit{Result} \\[4pt]
Dimensions match ($m=p$ and $n=q$) &
$c_{ij} = a_{ij} - b_{ij}$ \\
$A$ is a column vector; rows match ($m=p$; $n=1$) &
$c_{ij} = a_{i} - b_{ij}$ \\
$B$ is a column vector; rows match ($m=p$; $q=1$) &
$c_{ij} = a_{ij} - b_{i}$ \\
$A$ is a row vector; columns match ($m=1$; $n=q$) &
$c_{ij} = a_{j} - b_{ij}$ \\
$B$ is a row vector; columns match ($m=p$; $q=1$) &
$c_{ij} = a_{ij} - b_{j}$ \\
$A$ is a column vector; $B$ is a row vector ($n=1$; $p=1$) &
$c_{ij} = a_{i} - b_{j}$ \\
$A$ is a row vector; $B$ is a column vector ($m=1$; $q=1$) &
$c_{ij} = a_{j} - b_{i}$ \\
$A$ is a scalar ($m=1$ and $n=1$) &
$c_{ij} = a - b_{ij}$ \\
$B$ is a scalar ($p=1$ and $q=1$) &
$c_{ij} = a_{ij} - b$ \\
\end{tabular}
\end{center}
%
If none of the above conditions are satisfied the result is undefined
and an error is flagged.
Note that this convention makes it unnecessary, in most cases, to use
diagonal matrices to perform transformations by means of ordinary
matrix multiplication: if $Y = XV$, where $V$ is diagonal, it is
computationally much more convenient to obtain $Y$ via the instruction
%
\begin{code}
matrix Y = X .* v
\end{code}
%
where \texttt{v} is a row vector containing the diagonal of $V$.
In \textit{column-wise concatenation} of an $m\times n$ matrix $A$ and
an $m\times p$ matrix $B$, the result is an $m\times (n+p)$ matrix.
That is,
%
\begin{code}
matrix C = A ~ B
\end{code}
%
produces $C = \left[ \begin{array}{cc} A & B \end{array} \right]$.
\textit{Row-wise concatenation} of an $m\times n$ matrix $A$ and
an $p\times n$ matrix $B$ produces an $(m+p) \times n$ matrix.
That is,
%
\begin{code}
matrix C = A | B
\end{code}
%
produces $C = \left[ \begin{array}{cc} A \\ B \end{array} \right]$.
\section{Matrix--scalar operators}
\label{matrix-scalar-op}
For matrix $A$ and scalar $k$, the operators shown in
Table~\ref{tab:matrix-scalar-ops} are available. (Addition and
subtraction were discussed in section~\ref{matrix-op} but we include
them in the table for completeness.) In addition, for square $A$ and
integer $k \geq 0$, \verb|B = A^k| produces a matrix $B$ which is $A$
raised to the power $k$.
\begin{table}[htbp]
\centering
\begin{tabular}{ll}
\textit{Expression} & \textit{Effect} \\[4pt]
\texttt{matrix B = A * k} & $b_{ij} = k a_{ij}$ \\
\texttt{matrix B = A / k} & $b_{ij} = a_{ij} / k$ \\
\texttt{matrix B = k / A} & $b_{ij} = k / a_{ij}$ \\
\texttt{matrix B = A + k} & $b_{ij} = a_{ij} + k$ \\
\texttt{matrix B = A - k} & $b_{ij} = a_{ij} - k$ \\
\texttt{matrix B = k - A} & $b_{ij} = k - a_{ij}$ \\
\texttt{matrix B = A \% k} & $b_{ij} = a_{ij} \mbox{ modulo } k$ \\
\end{tabular}
\caption{Matrix--scalar operators}
\label{tab:matrix-scalar-ops}
\end{table}
\section{Matrix functions}
\label{sec:matrix-func}
Most of the \app{gretl} functions available for scalars and series
also apply to matrices in an element-by-element fashion, and as such
their behavior should be pretty obvious. This is the case for
functions such as \texttt{log}, \texttt{exp}, \texttt{sin}, etc.
These functions have the effects documented in relation to the
\texttt{genr} command. For example, if a matrix \texttt{A} is already
defined, then
%
\begin{code}
matrix B = sqrt(A)
\end{code}
%
generates a matrix such that $b_{ij} = \sqrt{a_{ij}}$. All such
functions require a single matrix as argument, or an expression which
evaluates to a single matrix.\footnote{Note that to find the ``matrix
square root'' you need the \texttt{cholesky} function (see below);
moreover, the \texttt{exp} function computes the exponential element
by element, and therefore does \emph{not} return the matrix
exponential unless the matrix is diagonal --- to get the matrix
exponential, use \texttt{mexp}.}
In this section, we review some aspects of \texttt{genr} functions that
apply specifically to matrices. A full account of each function is
available in the \GCR.
\newlength{\cwid}
\setlength{\cwid}{0.1\textwidth}
\begin{table}[htbp]
\centering
\input matfuncs.tex
\caption{Matrix functions by category}
\label{tab:matrix_funcs_cat}
\end{table}
\subsection{Matrix reshaping}
\label{matrix-mshape}
In addition to the methods discussed in sections
\ref{sec:matrix-create} and \ref{sec:matrix-sub}, a matrix can also be
created by re-arranging the elements of a pre-existing matrix. This is
accomplished via the \texttt{mshape} function. It takes three
arguments: the input matrix, $A$, and the rows and columns of the
target matrix, $r$ and $c$ respectively. Elements are read from $A$
and written to the target in column-major order. If $A$ contains
fewer elements than $n = r \times c$, they are repeated cyclically; if
$A$ has more elements, only the first $n$ are used.
For example:
\begin{code}
matrix a = mnormal(2,3)
a
matrix b = mshape(a,3,1)
b
matrix b = mshape(a,5,2)
b
\end{code}
produces
\begin{code}
? a
a
1.2323 0.99714 -0.39078
0.54363 0.43928 -0.48467
? matrix b = mshape(a,3,1)
Generated matrix b
? b
b
1.2323
0.54363
0.99714
? matrix b = mshape(a,5,2)
Replaced matrix b
? b
b
1.2323 -0.48467
0.54363 1.2323
0.99714 0.54363
0.43928 0.99714
-0.39078 0.43928
\end{code}
\subsection{Complex multiplication and division}
\label{sec:complex}
\app{Gretl} has no native provision for complex numbers. However,
basic operations can be performed on vectors of complex numbers by
using the convention that a vector of $n$ complex numbers is
represented as a $n \times 2$ matrix, where the first column contains
the real part and the second the imaginary part.
Addition and subtraction are trivial; the functions \texttt{cmult}
and \texttt{cdiv} compute the complex product and division,
respectively, of two input matrices, $A$ and $B$, representing complex
numbers. These matrices must have the same number of rows, $n$, and
either one or two columns. The first column contains the real part
and the second (if present) the imaginary part. The return value is
an $n \times 2$ matrix, or, if the result has no imaginary part, an
$n$-vector.
For example, suppose you have $z_1 = [ 1 + 2i , 3 + 4i ]'$ and $z_2 =
[ 1, i ]'$:
\begin{code}
? z1 = {1,2;3,4}
z1 = {1,2;3,4}
Generated matrix z1
? z2 = I(2)
z2 = I(2)
Generated matrix z2
? conj_z1 = z1 .* {1,-1}
conj_z1 = z1 .* {1,-1}
Generated matrix conj_z1
? eval cmult(z1,z2)
eval cmult(z1,z2)
1 2
-4 3
? eval cmult(z1,conj_z1)
eval cmult(z1,conj_z1)
5
25
\end{code}
\subsection{Multiple returns and the \texttt{null} keyword}
\label{matrix-multiples}
Some functions take one or more matrices as arguments and compute one
or more matrices; these are:
\begin{center}
\begin{tabular}{ll}
\texttt{eigensym} & Eigen-analysis of symmetric matrix \\
\texttt{eigengen} & Eigen-analysis of general matrix \\
\texttt{mols} & Matrix OLS \\
\texttt{qrdecomp} & QR decomposition \\
\texttt{svd} & Singular value decomposition (SVD)
\end{tabular}
\end{center}
The general rule is: the ``main'' result of the function is always
returned as the result proper. Auxiliary returns, if needed, are
retrieved using pre-existing matrices, which are passed to the
function as pointers (see \ref{funscope}). If such values are not
needed, the pointer may be substituted with the keyword \texttt{null}.
The syntax for \texttt{qrdecomp}, \texttt{eigensym} and
\texttt{eigengen} is of the form
%
\begin{code}
matrix B = func(A, &C)
\end{code}
%
The first argument, \texttt{A}, represents the input data, that is,
the matrix whose decomposition or analysis is required. The second
argument must be either the name of an existing matrix preceded by
\verb+&+ (to indicate the ``address'' of the matrix in question), in
which case an auxiliary result is written to that matrix, or the
keyword \texttt{null}, in which case the auxiliary result is not
produced, or is discarded.
In case a non-null second argument is given, the specified matrix will
be over-written with the auxiliary result. (It is not required that
the existing matrix be of the right dimensions to receive the result.)
The function \texttt{eigensym} computes the eigenvalues, and
optionally the right eigenvectors, of a symmetric $n \times n$ matrix.
The eigenvalues are returned directly in a column vector of length
$n$; if the eigenvectors are required, they are returned in an $n
\times n$ matrix. For example:
%
\begin{code}
matrix V
matrix E = eigensym(M, &V)
matrix E = eigensym(M, null)
\end{code}
%
In the first case \texttt{E} holds the eigenvalues of \texttt{M} and
\texttt{V} holds the eigenvectors. In the second, \texttt{E} holds
the eigenvalues but the eigenvectors are not computed.
The function \texttt{eigengen} computes the eigenvalues, and
optionally the eigenvectors, of a general $n \times n$ matrix. The
eigenvalues are returned directly in an $n \times 2$ matrix, the first
column holding the real components and the second column the imaginary
components.
If the eigenvectors are required (that is, if the second argument to
\texttt{eigengen} is not \texttt{null}), they are returned in an $n
\times n$ matrix. The column arrangement of this matrix is somewhat
non-trivial: the eigenvectors are stored in the same order as the
eigenvalues, but the real eigenvectors occupy one column, whereas
complex eigenvectors take two (the real part comes first); the total
number of columns is still $n$, because the conjugate eigenvector is
skipped. Example \ref{cmplx-evecs} provides a (hopefully) clarifying
example (see also subsection \ref{sec:complex}).
\begin{script}[htbp]
\caption{Complex eigenvalues and eigenvectors}
\label{cmplx-evecs}
\begin{scode}
set seed 34756
matrix v
A = mnormal(3,3)
/* do the eigen-analysis */
l = eigengen(A,&v)
/* eigenvalue 1 is real, 2 and 3 are complex conjugates */
print l
print v
/*
column 1 contains the first eigenvector (real)
*/
B = A*v[,1]
c = l[1,1] * v[,1]
/* B should equal c */
print B
print c
/*
columns 2:3 contain the real and imaginary parts
of eigenvector 2
*/
B = A*v[,2:3]
c = cmult(ones(3,1)*(l[2,]),v[,2:3])
/* B should equal c */
print B
print c
\end{scode}
\end{script}
The \texttt{qrdecomp} function computes the QR decomposition of an $m
\times n$ matrix $A$: $A = QR$, where $Q$ is an $m \times n$
orthogonal matrix and $R$ is an $n \times n$ upper triangular matrix.
The matrix $Q$ is returned directly, while $R$ can be retrieved via
the second argument. Here are two examples:
%
\begin{code}
matrix R
matrix Q = qrdecomp(M, &R)
matrix Q = qrdecomp(M, null)
\end{code}
%
In the first example, the triangular $R$ is saved as \texttt{R}; in
the second, $R$ is discarded. The first line above shows an example
of a ``simple declaration'' of a matrix: \texttt{R} is
declared to be a matrix variable but is not given any explicit value.
In this case the variable is initialized as a $1\times 1$ matrix whose
single element equals zero.
The syntax for \texttt{svd} is
%
\begin{code}
matrix B = func(A, &C, &D)
\end{code}
%
The function \texttt{svd} computes all or part of the singular value
decomposition of the real $m \times n$ matrix $A$. Let $k =
\mbox{min}(m, n)$. The decomposition is
\[
A = U \Sigma V'
\]
where $U$ is an $m \times k$ orthogonal matrix, $\Sigma$ is an $k
\times k$ diagonal matrix, and $V$ is an $k \times n$ orthogonal
matrix.\footnote{This is not the only definition of the SVD: some
writers define $U$ as $m \times m$, $\Sigma$ as $m \times n$ (with
$k$ non-zero diagonal elements) and $V$ as $n \times n$.} The
diagonal elements of $\Sigma$ are the singular values of $A$; they are
real and non-negative, and are returned in descending order. The
first $k$ columns of $U$ and $V$ are the left and right singular
vectors of $A$.
The \texttt{svd} function returns the singular values, in a vector of
length $k$. The left and/or right singular vectors may be obtained by
supplying non-null values for the second and/or third arguments
respectively. For example:
%
\begin{code}
matrix s = svd(A, &U, &V)
matrix s = svd(A, null, null)
matrix s = svd(A, null, &V)
\end{code}
%
In the first case both sets of singular vectors are obtained, in the
second case only the singular values are obtained; and in the third,
the right singular vectors are obtained but $U$ is not computed.
\emph{Please note}: when the third argument is non-null, it is
actually $V'$ that is provided. To reconstitute the original matrix
from its SVD, one can do:
%
\begin{code}
matrix s = svd(A, &U, &V)
matrix B = (U.*s)*V
\end{code}
%
Finally, the syntax for \texttt{mols} is
%
\begin{code}
matrix B = mols(Y, X, &U)
\end{code}
%
This function returns the OLS estimates obtained by regressing the $T
\times n$ matrix $Y$ on the $T \times k$ matrix $X$, that is, a $k
\times n$ matrix holding $(X'X)^{-1} X'Y$. The Cholesky decomposition
is used. The matrix $U$, if not \texttt{null}, is used to store the
residuals.
\subsection{Reading and writing matrices from/to text files}
\label{sec:matrix-csv}
The two functions \texttt{mread} and \texttt{mwrite} can be used for
basic matrix input/output. This can be useful to enable \app{gretl} to
exchange data with other programs.
The \texttt{mread} function accepts one string parameter: the name of
the (plain text) file from which the matrix is to be read. The file
in question must conform to the following rules:
%
\begin{enumerate}
\item The columns must be separated by spaces or tab characters.
\item The decimal separator must be the dot ``\texttt{.}'' character.
\item The first line in the file must contain two integers, separated
by a space or a tab, indicating the number of rows and columns,
respectively.
\end{enumerate}
Should an error occur (such as the file being badly formatted or
inaccessible), an empty matrix (see section~\ref{sec:emptymatrix}) is
returned.
The complementary function \texttt{mwrite} produces text files
formatted as described above. The column separator is the tab
character, so import into spreadsheets should be straightforward.
Usage is illustrated in example \ref{matrix-rw}. Matrices stored via
the \texttt{mwrite} command can be easily read by other programs; the
following table summarizes the appropriate commands for reading a
matrix \texttt{A} from a file called \texttt{a.mat} in some
widely-used programs.\footnote{Matlab users may find the Octave
example helpful, since the two programs are mostly compatible with one
another.}
\begin{center}
\begin{tabular}{rl}
\textbf{Program} & \textbf{Sample code} \\
\hline
GAUSS & \verb|tmp[] = load a.mat;| \\
& \verb|A = reshape(tmp[3:rows(tmp)],tmp[1],tmp[2]);| \\
Octave & \verb|fd = fopen("a.mat");| \\
& \verb|[r,c] = fscanf(fd, "%d %d", "C");| \\
& \verb|A = reshape(fscanf(fd, "%g", r*c),c,r)';| \\
& \verb|fclose(fd);| \\
Ox & \verb|decl A = loadmat("a.mat");| \\
R & \verb|A <- as.matrix(read.table("a.mat", skip=1))| \\
\hline
\end{tabular}
\end{center}
\begin{script}[htbp]
\caption{Matrix input/output via text files}
\label{matrix-rw}
\begin{scode}
nulldata 64
scalar n = 3
string f1 = "a.csv"
string f2 = "b.csv"
matrix a = mnormal(n,n)
matrix b = inv(a)
err = mwrite(a, f1)
if err != 0
fprintf "Failed to write %s\n", f1
else
err = mwrite(b, f2)
endif
if err != 0
fprintf "Failed to write %s\n", f2
else
c = mread(f1)
d = mread(f2)
a = c*d
printf "The following matrix should be an identity matrix\n"
print a
endif
\end{scode}
\end{script}
\section{Matrix accessors}
\label{matrix-accessors}
In addition to the matrix functions discussed above,
various ``accessor'' strings allow you to create copies of internal
matrices associated with models previously estimated.
These are set out in Table~\ref{tab:matrix-accessors}.
\begin{table}[htbp]
\centering
\begin{tabular}{ll}
\texttt{\$coeff} & vector of estimated coefficients \\
\texttt{\$compan} & companion matrix (after VAR or VECM estimation) \\
\texttt{\$jalpha} & matrix $\alpha$ (loadings) from Johansen's procedure \\
\texttt{\$jbeta} & matrix $\beta$ (cointegration vectors) from
Johansen's procedure \\
\texttt{\$jvbeta} & covariance matrix for the unrestricted elements of
$\beta$ from Johansen's procedure \\
\texttt{\$rho} & autoregressive coefficients for error process \\
\texttt{\$sigma} & residual covariance matrix \\
\texttt{\$stderr} & vector of estimated standard errors \\
\texttt{\$uhat} & matrix of residuals \\
\texttt{\$vcv} & covariance matrix of parameter estimates \\
\texttt{\$yhat} & matrix of fitted values
\end{tabular}
\caption{Matrix accessors for model data}
\label{tab:matrix-accessors}
\end{table}
Many of the accessors in Table~\ref{tab:matrix-accessors} behave
somewhat differently depending on the sort of model that is
referenced, as follows:
\begin{itemize}
\item Single-equation models: \texttt{\$sigma} gets a scalar (the
standard error of the residuals); \texttt{\$uhat} and
\texttt{\$yhat} get series.
\item All system estimators: \texttt{\$sigma} gets the cross-equation
residual covariance matrix, \texttt{\$uhat} gets a matrix of
residuals, one column per equation.
\item VARs and VECMs: \texttt{\$stderr} and \texttt{\$yhat} are not
available; \texttt{\$coeff} gets a matrix of coefficients, one
column per equation.
\end{itemize}
If the accessors are given without any prefix, they retrieve results
from the last model estimated, if any. Alternatively, they may be
prefixed with the name of a saved model plus a period (\texttt{.}), in
which case they retrieve results from the specified model. Here are
some examples:
%
\begin{code}
matrix u = $uhat
matrix b = m1.$coeff
matrix v2 = m1.$vcv[1:2,1:2]
\end{code}
%$
The first command grabs the residuals from the last model; the second
grabs the coefficient vector from model \texttt{m1}; and the third
(which uses the mechanism of sub-matrix selection described above)
grabs a portion of the covariance matrix from model \texttt{m1}.
If the model in question a VAR or VECM (only) \verb|$compan| returns
the companion matrix.
After a vector error correction model is estimated via Johansen's
procedure, the matrices \verb|$jalpha| and \verb|$jbeta| are
also available. These have a number of columns equal to the chosen
cointegration rank; therefore, the product
\begin{code}
matrix Pi = $jalpha * $jbeta'
\end{code}
returns the reduced-rank estimate of $A(1)$. Since $\beta$ is
automatically identified via the Phillips normalization (see section
\ref{sec:johansen-ident}), its unrestricted elements do have a proper
covariance matrix, which can be retrieved through the
\verb|$jvbeta| accessor.
\section{Namespace issues}
\label{matrix-namespace}
Matrices share a common namespace with data series and scalar
variables. In other words, no two objects of any of these types can
have the same name. It is an error to attempt to change the type of
an existing variable, for example:
%
\begin{code}
scalar x = 3
matrix x = ones(2,2) # wrong!
\end{code}
%
It is possible, however, to delete or rename an existing variable then
reuse the name for a variable of a different type:
\begin{code}
scalar x = 3
delete x
matrix x = ones(2,2) # OK
\end{code}
\section{Creating a data series from a matrix}
\label{matrix-create-series}
Section~\ref{sec:matrix-create} above describes how to create a matrix
from a data series or set of series. You may sometimes wish to go in
the opposite direction, that is, to copy values from a matrix
into a regular data series. The syntax for this operation is
%
\begin{textcode}
series \textsl{sname} = \textsl{mspec}
\end{textcode}
%
where \ttsl{sname} is the name of the series to create and
\ttsl{mspec} is the name of the matrix to copy from, possibly followed
by a matrix selection expression. Here are two examples.
%
\begin{code}
series s = x
series u1 = U[,1]
\end{code}
%
It is assumed that \texttt{x} and \texttt{U} are pre-existing
matrices. In the second example the series \texttt{u1} is formed from
the first column of the matrix \texttt{U}.
For this operation to work, the matrix (or matrix selection) must be a
vector with length equal to either the full length of the current
dataset, $n$, or the length of the current sample range, $n^{\prime}$.
If $n^{\prime} < n$ then only $n^{\prime}$ elements are drawn from the
matrix; if the matrix or selection comprises $n$ elements, the
$n^{\prime}$ values starting at element $t_1$ are used, where $t_1$
represents the starting observation of the sample range. Any values
in the series that are not assigned from the matrix are set to the
missing code.
\section{Matrices and lists}
\label{matrix-and-list}
To facilitate the manipulation of named lists of variables (see
Chapter~\ref{chap-persist}), it is possible to convert between
matrices and lists. In section~\ref{sec:matrix-create} above we mentioned
the facility for creating a matrix from a list of variables, as in
%
\begin{code}
matrix M = { listname }
\end{code}
%
That formulation, with the name of the list enclosed in braces, builds
a matrix whose columns hold the variables referenced in the list.
What we are now describing is a different matter: if we say
%
\begin{code}
matrix M = listname
\end{code}
%
(without the braces), we get a row vector whose elements are
the ID numbers of the variables in the list. This special case
of matrix generation cannot be embedded in a compound
expression. The syntax must be as shown above, namely simple
assignment of a list to a matrix.
To go in the other direction, you can include a matrix on the
right-hand side of an expression that defines a list, as in
%
\begin{code}
list Xl = M
\end{code}
%
where \texttt{M} is a matrix. The matrix must be suitable for
conversion; that is, it must be a row or column vector containing
non-negative whole-number values, none of which exceeds the highest ID
number of a variable (series or scalar) in the current dataset.
Example~\ref{normalize-list} illustrates the use of this sort of
conversion to ``normalize'' a list, moving the constant (variable 0)
to first position.
\begin{script}[htbp]
\caption{Manipulating a list}
\label{normalize-list}
\begin{scode}
function normalize_list (matrix *x)
# If the matrix (representing a list) contains var 0,
# but not in first position, move it to first position
if (x[1] != 0)
scalar k = cols(x)
loop for (i=2; i<=k; i++) --quiet
if (x[i] = 0)
x[i] = x[1]
x[1] = 0
break
endif
end loop
end if
end function
open data9-7
list Xl = 2 3 0 4
matrix x = Xl
normalize_list(&x)
list Xl = x
\end{scode}
\end{script}
\section{Deleting a matrix}
\label{matrix-delete}
To delete a matrix, just write
%
\begin{code}
delete M
\end{code}
%
where \texttt{M} is the name of the matrix to be deleted.
\section{Printing a matrix}
To print a matrix, the easiest way is to give the name of the matrix
in question on a line by itself, which is equivalent to using the
\cmd{print} command:
%
\begin{code}
matrix M = mnormal(100,2)
M
print M
\end{code}
You can get finer control on the formatting of output by using the
\cmd{printf} command: for example, the following code
%
\begin{code}
matrix Id = I(2)
printf "%10.3f", Id
\end{code}
%
produces
\begin{code}
? print Id
print Id
Id (2 x 2)
1 0
0 1
? printf "%10.3f", Id
1.000 0.000
0.000 1.000
\end{code}
For presentation purposes you may wish to give titles to the columns
of a matrix. For this you can use the \cmd{colnames} function: the first
argument is a matrix and the second is either a named list of variables,
whose names will be used as headings, or a string that contains as many
space-separated substrings as the matrix has columns. For example,
%
\begin{code}
? matrix M = mnormal(3,3)
? colnames(M, "foo bar baz")
? print M
M (3 x 3)
foo bar baz
1.7102 -0.76072 0.089406
-0.99780 -1.9003 -0.25123
-0.91762 -0.39237 -1.6114
\end{code}
\section{Example: OLS using matrices}
\label{matrix-example}
Example \ref{matrixOLS} shows how matrix methods can be used to
replicate gretl's built-in OLS functionality.
\begin{script}[htbp]
\caption{OLS via matrix methods}
\label{matrixOLS}
\begin{scode}
open data4-1
matrix X = { const, sqft }
matrix y = { price }
matrix b = invpd(X'X) * X'y
print "estimated coefficient vector"
b
matrix u = y - X*b
scalar SSR = u'u
scalar s2 = SSR / (rows(X) - rows(b))
matrix V = s2 * inv(X'X)
V
matrix se = sqrt(diag(V))
print "estimated standard errors"
se
# compare with built-in function
ols price const sqft --vcv
\end{scode}
\end{script}
