\chapter{Matrix manipulation}
\label{chap:matrices}

Together with the other two basic types of data (series and scalars),
gretl offers a quite comprehensive array of matrix methods. This
chapter illustrates the peculiarities of matrix syntax and discusses
briefly some of the more advanced matrix functions. For a full listing
of matrix functions and a comprehensive account of their syntax,
please refer to the \GCR.

In this chapter we're concerned with real matrices; most of the points
made here also apply to complex matrices but see the following chapter
for additional specifics on the complex case.

\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---either in numerical form, or 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 Via a suitable expression that references 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 series are separated by commas.  Besides names
of existing series, you can use expressions that yield a series
result.  For example, given a series \texttt{x} you could do
%
\begin{code}
matrix A = {x, x^2}
\end{code}
%
Each series occupies a column (and there can only be one series per
column).  The semicolon cannot be used 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.

Instead of giving an explicit list of series, you may instead provide
the \textit{name of a saved list} (see
Chapter~\ref{chap:lists-strings}), 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 series,
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 (series 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 series using horizontal concatenation (see below), as in
%
\begin{code}
matrix A = {const}~{dataset}
\end{code}
%

By default, when you build a matrix from series that include missing
values any observations (rows) that contain \texttt{NA}s are skipped.
This behavior can be modified via the command \dtk{set skip_missing
  off}.  In that case missing values are represented in the matrix as
NaN (``Not a Number'').  Note that as per the IEEE standard,
arithmetic operations involving one or more NaNs always produce a NaN
result. Alternatively, you can take greater control over the
observations (data rows) that are included in the matrix using the
``set'' variable \dtk{matrix_mask}, as in
%
\begin{code}
set matrix_mask msk
\end{code}
%
where \texttt{msk} is the name of a series.  Subsequent commands that
form matrices from series or lists will include only observations
for which \texttt{msk} has non-zero (and non-missing) values. You
can remove this mask via the command \dtk{set matrix_mask null}.

\tip{Names of matrices must satisfy the same requirements as names of
  gretl variables in general: the name can be no longer than 31
  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 syntax 
%
\begin{code}
matrix A = {}
\end{code}
%
creates an empty matrix---a matrix with zero rows and zero columns.

The main purpose of this construct is to allow 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}
%
produce a matrix \texttt{B} that is identical to \texttt{X}.

\begin{table}[htbp]
\centering
\begin{tabular}{lc@{\hspace{5em}}lc}
\textit{Function} & \textit{Return value} & \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} \\
  \texttt{diag(A)} & \verb|{}| &
  \texttt{vec(A)} & \verb|{}| \\
  \texttt{vech(A)} & \verb|{}| &
  \texttt{unvech(A)} & \verb|{}| \\
\end{tabular}
\caption{Valid functions on an empty matrix, \texttt{A}}
\label{tab:empty-matrix-funcs}
\end{table}

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 function calls with an empty matrix argument are listed in
Table~\ref{tab:empty-matrix-funcs}; other matrix functions generate an
error if given an empty argument.  In line with the above
interpretation, the functions \texttt{I}, \texttt{ones},
\texttt{zeros}, \texttt{mnormal} and \texttt{muniform} return an empty
matrix when or more of the arguments is 0, as well as the function
\texttt{nullspace} when its argument is of full column rank.

While the $0 \times 0$ matrix produced by \verb|{}| is the most common
sort of empty matrix, matrices with one zero and one positive
dimension are also supported. For example, the call
\begin{code}
  matrix A = zeros(3,0)
\end{code}
produces a $3 \times 0$ matrix, as does
\begin{code}
matrix A = zeros(3,1)
A = A[,-1] # delete column 1 (see below)
\end{code}
The handling of such matrices by gretl is very similar to that in
\textsf{Matlab}.

\section{Selecting submatrices}
\label{sec:matrix-sub}

You can select submatrices 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 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}.

If the matrix in question is $n\times 1$ or $1\times m$, it is
OK to give just one index specifier and omit the comma. For example,
\texttt{A[2]} selects the second element of \texttt{A} if \texttt{A}
is a vector. Otherwise the comma is mandatory.

In addition there are some predefined index specifications,
represented by the keywords \texttt{diag}, \texttt{lower},
\texttt{upper}, \texttt{real}, \texttt{imag} and \texttt{end}. With
the exception of \texttt{end}, these keywords imply specific row
\textit{and} column selections, and therefore cannot be combined with
any additional, comma-separated term.

\begin{itemize}
\item The \texttt{diag} specification selects the principal diagonal
  of a matrix.
\item \texttt{lower} and \texttt{upper} select, respectively, the
  elements of a matrix below and those above the principal diagonal.
\item \texttt{real} and \texttt{imag} are specific to complex matrices
  and are described in chapter~\ref{chap:complex}.
\item \texttt{end} selects the last element in a given row or column.
  It can be employed in arithmetical expressions, so for example
  \texttt{end-1} accesses the second-last element in a row or column.
\end{itemize}

Submatrix selections cane be used 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$, then the lower triangle of $A$:
%
\begin{code}
matrix A = {1, 2, 3; 4, 5, 6; 7, 8, 9}
matrix B = A[1:2,2:3]
matrix C = A[lower]
\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}

When the \texttt{lower} and \texttt{upper} selections are used on the
right, they yield a vector holding the elements in their scope. The
ordering of the elements is column-major in both cases, as illustrated
below for the $4 \times 4$ case.
\begin{center}
  \begin{tabular}{cccc}
    $d$ & 1 & 2 & 4 \\
    1 & $d$ & 3 & 5 \\
    2 & 4 & $d$ & 6 \\
    3 & 5 & 6 & $d$
  \end{tabular}
\end{center}
This means that \texttt{lower} and \texttt{upper} do not produce the
same result for symmetric matrices bigger than $3 \times 3$, which may
seem unfortunate, but it gives the user a degree of flexibility in
respect of the ordering of the elements. Suppose you have a
non-symmetric matrix $M$ and you'd like to extract the infradiagonal
elements in \textit{row}-major order: \texttt{(M')[upper]} will do the
job.

When \texttt{lower} and \texttt{upper} are used on the left, the
replacement must be either (a) a vector of length equal to the number
of elements in the selection or (b) a scalar value. In case (a) the
elements of the target matrix are filled in column-major order; in
case (b) they are all set using the scalar.

One possible use of these tools is taking (say) a lower triangular
matrix and rendering it symmetric by copying the elements from beneath
the diagonal to above. The way to get this right (assuming you have a
lower triangular matrix \texttt{L}) is
\begin{code}
L[upper] = (L')[upper]  # note: not L[upper] = L[lower]
\end{code}

\section{Deleting rows or columns}
\label{sec:neg-indices}

A variant of submatrix notation is available for convenience in
dropping specified rows and/or columns from a matrix, namely giving
negative values for the indices. Here is a simple example,
%
\begin{code}
matrix A = {1, 2, 3; 4, 5, 6; 7, 8, 9}
matrix B = A[-2,-3]
\end{code}
%
which creates \texttt{B} as a $2\times 2$ matrix which drops row 2 and
column 3 from \texttt{A}. Negative indices can also be given in the
form of an index vector:
%
\begin{code}
matrix rdrop = {-1,-3,-5}
matrix B = A[rdrop,]
\end{code}
%
In this case \texttt{B} is formed by dropping rows 1, 3 and 5 from
\texttt{A} (which must have at least 5 rows), but retaining the column
dimension of \texttt{A}.

There are two limitations on the use of negative indices. First, the
\texttt{from:to} range syntax described in the previous section is not
available, but you can use the \texttt{seq} function to achieve an
equivalent effect, as in
%
\begin{code}
matrix A = muniform(1, 10)
matrix B = A[,-seq(3,7)]
\end{code}
%
where \texttt{B} drops columns 3 to 7 from \texttt{A}. Second, use of
negative indices is valid only on the right-hand side of a matrix
calculation; there is no ``negative index'' equivalent of assignment
to a submatrix, as in
%
\begin{code}
A[1:3,] = ones(3, cols(A))
\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 \\
\verb|\|    & matrix ``left division'' (see below) \\
\texttt{/}  & matrix ``right division'' (see below) \\
\verb+~+    & column-wise concatenation \\
\verb+|+    & row-wise concatenation \\
\texttt{**} & Kronecker product \\
\texttt{==}  & test for equality \\
\texttt{!=} & test for inequality
\end{tabular}
\end{center}

In addition, the following operators (``dot'' operators) apply on an
element-by-element basis:

\begin{center}
\begin{tabular}{ccccccccccc}
\texttt{.+}  &  \texttt{.-}  &
\texttt{.*}  &  \texttt{./}  &  \verb+.^+  &
\texttt{.=}  &  \texttt{.>}  &  \texttt{.<} &
\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 syntax. In the special case $X = Y$, however, the two are not
exactly equivalent. The former expression uses a specialized algorithm
with two advantages: it is more efficient computationally, and ensures
that the result is free of machine precision artifacts that may render
it numerically non-symmetric. This, however, is unlikely to be an
issue unless your $X$ matrix is rather large (at least several
hundreds rows/columns).

In matrix ``left division'', the statement 
%
\begin{code}
matrix X = A \ B
\end{code}
%
is interpreted as a request to find the matrix $X$ that solves $AX=B$,
so $A$ and $B$ must have the same number of rows. If $A$ is a square
matrix, this is in principle equivalent to $A^{-1}B$, which fails if
$A$ is singular; the numerical method employed here is the LU
decomposition.  If $A$ is a $T \times k$ matrix with $T > k$, then $X$
is the least-squares solution, $X = (A'A)^{-1}A'B$, which fails if
$A'A$ is singular; the numerical method is the QR decomposition.
Otherwise, the operation fails.

For matrix ``right division'', as in \texttt{X = A / B}, $X$ is the
matrix that solves $XB = A$, so $A$ and $B$ must have the same number
of columns. If $B$ is non-singular this is in principle equivalent to
$AB^{-1}$, otherwise $X$ is the least-squares solution.

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$.

\subsection{Concatenation}

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]$.

In general the two matrix operands must be conformable, as indicated
above. But as a special case for convenience, a scalar or $1 \times 1$
matrix can be concatenated with a matrix of any size; the single value
is automatically repeated, if necessary, to form a row or column
vector of length conformable with the matrix in question. So for
example
%
\begin{code}
matrix A = {1, 2, 3}
matrix B = A | 1
\end{code}
%
produces $B = \left[ \begin{array}{ccc} 1 & 2 & 3 \\
  1 & 1 & 1 \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
scalar $x$, \verb|B = A^x| produces a matrix $B$ which is $A$ raised
to the power $x$, but only if either of two conditions are
satisfied. First, if $x$ is a non-negative integer then Golub and Van
Loan's ``Binary Powering'' Algorithm 11.2.2 is used---see
\cite{golub96}---and $A$ can then be a generic square matrix. Second,
if $A$ is positive semidefinite the power is computed via its
eigen-decomposition and $x$ can be a real number, subject to the
constraint that $x$ can be negative only if $A$ is invertible.

\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 functions available for scalars and series also apply to
matrices on an element-by-element basis. This is the case for
\texttt{log}, \texttt{exp}, \texttt{sqrt}, \texttt{sin} and many
others. 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).
  And since the \texttt{exp} function computes the exponential element
  by element, it 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 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{Multiple returns and the \texttt{null} keyword}
\label{matrix-multiples}

The functions listed below take one or more matrices as arguments and
compute one or more matrices.

\begin{center}
\begin{tabular}{ll}
\texttt{eigensym} & Eigen-analysis of symmetric matrix \\
\texttt{eigen}    & Eigen-analysis of general matrix \\
\texttt{mols}     & Matrix OLS \\
\texttt{qrdecomp} & QR decomposition \\
\texttt{svd}      & Singular value decomposition (SVD) 
\end{tabular}
\end{center}

The ``main'' result of the function is always returned as the result
proper. Auxiliary returns, if wanted, are retrieved using pre-existing
matrices, passed to the function in ``pointer'' form
(Section~\ref{sec:funscope}). If such values are not needed, the
pointer may be substituted with the keyword \texttt{null}.

The syntax for \texttt{qrdecomp} and \texttt{eigensym} 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.

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{eigen} computes the eigenvalues, and optionally
the right and/or left eigenvectors, of a general $n \times n$
matrix. Following the input matrix argument there are two slots for
matrix-addresses, the first to retrieve the right eigenvectors and the
second for the left.  Calls to this function should therefore conform
to one of the following patterns.
\begin{code}
# get the eigenvalues only
matrix E = eigen(M)

# get the right eigenvectors as well
matrix V = {}
matrix E = eigen(M, &V)

# get both sets of eigenvectors
matrix V = {}
matrix W = {}
matrix E = eigen(M, &V, &W)

# get the left eigenvectors but not the right
matrix W = {}
matrix E = eigen(M, null, &W)
\end{code}

The eigenvalues are returned directly in a complex $n$-vector. If the
eigenvectors are wanted they are returned in a $n \times n$ complex
matrix.

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 but is not given any explicit value. The result is that it is
automatically initialized as an empty matrix (see
Section~\ref{sec:emptymatrix}).

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 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 may start with any number of comment lines, defined as
lines that start with the hash mark, ``\texttt{\#}''; such lines are
ignored.  Beyond that, the content must conform to the following
rules:
%
\begin{enumerate}
\item The first non-comment line must contain two integers, separated
  by a space or a tab, indicating the number of rows and columns,
  respectively.
\item The columns must be separated by spaces or tab characters.
\item The decimal separator must be the dot ``\texttt{.}'' character.
\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 Listing~\ref{ex:matrix-files}.  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.} Note that the Python example requires that the
\texttt{numpy} module is loaded.

\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))| \\
    Python & \verb|A = numpy.loadtxt('a.mat', skiprows=1)| \\
    Julia  & \verb|A = readdlm("a.mat", skipstart=1)| \\
  \hline
\end{tabular}
\end{center}


\begin{script}[htbp]
  \scriptinfo{matrix-files}{Matrix input/output via text files}
  \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}

Optionally, the \texttt{mwrite} and \texttt{mread} functions can use
gzip compression: this is invoked if the name of the matrix file has
the suffix ``\texttt{.gz}.'' In this case the elements of the matrix
are written in a single column. Note, however, that compression should
not be applied when writing matrices for reading by third-party
software unless you are sure that the software can handle compressed
data.

\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}
  \dollar{coeff}  & matrix of estimated coefficients \\
  \dollar{compan} & companion matrix (after VAR or VECM estimation) \\
  \dollar{jalpha} & matrix $\alpha$ (loadings) from Johansen's procedure \\
  \dollar{jbeta}  & matrix $\beta$ (cointegration vectors) from
  Johansen's procedure \\
  \dollar{jvbeta} & covariance matrix for the unrestricted elements of 
  $\beta$ from Johansen's procedure \\
  \dollar{rho}    & autoregressive coefficients for error process \\
  \dollar{sigma}  & residual covariance matrix \\
  \dollar{stderr} & matrix of estimated standard errors \\
  \dollar{uhat}   & matrix of residuals \\
  \dollar{vcv}    & covariance matrix of parameter estimates \\
  \dollar{vma}    & VMA matrices in stacked form (see section
    \ref{sec:var-estim}) \\
  \dollar{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: \dollar{sigma} gets a scalar (the
  standard error of the regression); \dollar{coeff} and
  \dollar{stderr} get column vectors; \dollar{uhat} and
  \dollar{yhat} get series.
\item System estimators: \dollar{sigma} gets the cross-equation
  residual covariance matrix; \dollar{uhat} and \dollar{yhat} get
  matrices with one column per equation.  The format of \dollar{coeff}
  and \dollar{stderr} depends on the nature of the system: for VARs
  and VECMs (where the matrix of regressors is the same for all
  equations) these return matrices with one column per equation, but
  for other system estimators they return a big column vector.
\item VARs and VECMs: \dollar{vcv} is not available, but 
  $X'X^{-1}$ (where $X$ is the common matrix of regressors) is
  available as \dollar{xtxinv}, 
  such that for VARs and VECMs (without restrictions on $\alpha$) a vcv
  equivalent can be easily and efficiently constructed as 
  \dollar{sigma} ** \dollar{xtxinv}.
\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 submatrix selection described above)
grabs a portion of the covariance matrix from model \texttt{m1}.

If the model in question a VAR or VECM (only) \dollar{compan} and
\dollar{vma} return the companion matrix and the VMA matrices in
stacked form, respectively (see section \ref{sec:var-estim} for
details).  After a vector error correction model is estimated via
Johansen's procedure, the matrices \dollar{jalpha} and \dollar{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
\dollar{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 series (see
Chapter~\ref{chap:lists-strings}), 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 series 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 series 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 integer values, none of which exceeds the highest ID
number of a series in the current dataset.

Listing~\ref{ex: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]
  \scriptinfo{normalize-list}{Manipulating a list}
\begin{scode}
function void 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++)
        if (x[i] == 0)
            x[i] = x[1]
            x[1] = 0
            break
         endif
     endloop
  endif
end function

open data9-7
list Xl = 2 3 0 4
matrix x = Xl
normalize_list(&x)
list Xl = x
list Xl print
\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, as illustrated in the interactive session below:
%
\begin{code}
? matrix Id = I(2)
 matrix Id = I(2)
Generated matrix Id
? 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 function \cmd{cnameset}. The
first argument to this function is a matrix and the second can take
any of three forms: a string that contains as many space-separated
substrings as the matrix has columns; an array of strings (see
Section~\ref{sec:arrays}) of the right length; or a named list of
series, whose names will be used as headings.  For example,
%
\begin{code}
? matrix M = mnormal(3,3)
? cnameset(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}
%
Row names may be added via the function \cmd{rnameset}, which works in
the same way.

\section{Example: OLS using matrices}
\label{matrix-example}

Listing~\ref{ex:matrix-ols} shows how matrix methods can be used to
replicate gretl's built-in OLS functionality.

\begin{script}[htbp]
  \scriptinfo{matrix-ols}{OLS via matrix methods}
\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}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "gretl-guide"
%%% End: 
