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