\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 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 semicolon 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 variables, you may instead
provide the \textit{name of a saved list} (see
Chapter~\ref{chap:persistent}), 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}
%

By default, when you build a matrix from series that include missing
values the data rows that contain \texttt{NA}s are skipped.  But you
can modify this behavior via the command \texttt{set skip\_missing
  off}.  In that case \texttt{NA}s are converted to NaN (``Not a
Number'').  In the IEEE floating-point standard, arithmetic operations
involving NaN always produce NaN. Alternatively, you can take greater
control over the observations (data rows) that are included in 
the matrix using the ``set'' variable \texttt{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 \texttt{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 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}.


\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} \\
\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 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.

\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 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 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{Deleting rows or columns}
\label{sec:neg-indices}

A variant of sub-matrix 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 sub-matrix, 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. Consider, however, that in the special case $X = Y$, the two
are not exactly equivalent: the former expression uses a specialized,
optimized algorithm which has the double advantage of being more
efficient computationally and of ensuring that the result will be free
by construction of machine precision artifacts that may render it
numerically non-symmetric. This, however, is unlikely to affect you
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$.
If $B$ 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 employed here is
the QR decomposition.  Otherwise, the operation necessarily fails.

For matrix ``right division'', as in \texttt{X = A / B}, $X$ is the
matrix that solves $XB = A$, in principle equivalent to $AB^{-1}$.

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 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}

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. Listing \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 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 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.} 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]
  \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}

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}.
\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) \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 variables (see
Chapter~\ref{chap:persistent}), 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 integer values, none of which exceeds the highest ID
number of a series in the current dataset.

Listing~\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 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++) --quiet
        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
\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 \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}

Listing~\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}


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