\SweaveOpts{results=hide, prefix.string=vigpics/sur}
\include{zinput}
%\VignetteIndexEntry{Seemingly Unrelated Regression}
%\VignetteDepends{Zelig}
%\VignetteKeyWords{model,continuous, least squares,error}
%\VignettePackage{Zelig}
\begin{document}
\nobibliography*
<<beforepkgs, echo=FALSE>>=
 before=search()
@
<<loadLibrary, echo=F,results=hide>>=
library(Zelig)
@ 

\section{{\tt sur}: Seemingly Unrelated Regression}
\label{sur}

\texttt{sur} extends ordinary least squares analysis to estimate system of linear 
equations with correlated error terms. The seemingly unrelated regression model can
be viewed as a special case of generalized least squares.
\subsubsection{Syntax}
\begin{verbatim}
> fml <- list ("mu1" = Y1 ~ X1,
               "mu2" = Y2 ~ X2,
               "mu3" = Y3 ~ X3)
> z.out<-zelig(formula = fml, model = "2sls", data = mydata)
> x.out <- setx(z.out)
> s.out <- sim(z.out, x = x.out)
\end{verbatim}
\subsubsection{Inputs}
\texttt{sur} regression specification has at least $M$ equations 
($M \ge 2$) corresponding to the dependent variables ($Y_1, Y_2, \ldots, Y_M$).
\begin{itemize}
\item \texttt{formula}:a list whose elements are formulae corresponding to 
the $M$ equations and their respective dependent and explanatory variables.
For example, when there are no constraints on the coefficients:
<<Inputs.list1>>=
 fml <- list ("mu1" = Y1 ~ X1,
               "mu2" = Y2 ~ X2,
               "mu3" = Y3 ~ X3)
@ 
\texttt{"mu1"} is the label for the first equation with Y1 as the dependent variable
and X1 as the explanatory variable. Similarly \texttt{"mu2"} and \texttt{"mu3"} are the
labels for the Y2 and Y3 equations.
\item \texttt{tag}: Users can also put constraints on the coefficients by using
the special function \texttt{tag}. \texttt{tag} takes two parameters. The first
parameter is the variable whose coefficient needs to be constrained and the second
parameter is label for the constrained coefficient. Each label uniquely identifies
the constrained coefficient. For example:
<<Inputs.list2>>=
 fml <- list ("mu1" = Y1 ~ tag(Xc,"constrain1")+ X1,
               "mu2" = Y2 ~ tag(Xc,"constrain1")+ X2,
               "mu3" = Y3 ~ X3)
@ 
\end{itemize}
\subsubsection{Additional Inputs}
\texttt{sur} takes the following additional inputs for model
specifications:
\begin{itemize}
\item \texttt{TX}: an optional matrix to transform the regressor
matrix and, hence, also the coefficient vector (see details). Default is \texttt{NULL}.
\item \texttt{maxiter}: maximum number of iterations.
\item \texttt{tol}: tolerance level indicating when to stop the iteration.
\item \texttt{rcovformula}: formula to calculate the estimated residual covariance
matrix (see details). Default is equal to 1.
\item \texttt{probdfsys}: use the degrees of freedom of the whole system
(in place of the degrees of freedom of the single equation to calculate probability
values for the t-test of individual parameters. 
\item \texttt{solvetol}: tolerance level for detecting linear dependencies when 
inverting a matrix or calculating a determinant. Default is \texttt {solvetol}= \begin{verbatim}.Machine\$double.eps.\end{verbatim}
\item \texttt{saveMemory}: logical. Save memory by omitting some calculation that are
not crucial for the basic estimate (e.g McElroy's $R^2$).
\end{itemize}
\subsubsection{Details}
The matrix \texttt{TX} transforms the regressor matrix ($X$) by $X\ast=X \times TX$. Thus,
the vector of coefficients is now $b=TX \times b\ast$ where $b$ is the original(stacked) 
vector of all coefficients and $b\ast$ is the new coefficient vector that is estimated instead.
Thus, the elements of vector $b$ and $b_i = \sum_j TX_{ij}\times b_j\ast$. The $TX$ matrix can be
used to change the order of the coefficients and also to restrict coefficients (if $TX$ has 
less columns than it has rows). 
If iterated (with \texttt{maxit}>1), the covergence criterion is
\begin{eqnarray*}
\sqrt{\frac{\sum_i(b_{i,g}-b_{i,g-1})^2}{\sum_ib_{i,g-1}^2}}< tol
\end{eqnarray*}
where $b_{i,g}$ is the ith coefficient of the gth iteration step.
The formula (\texttt{rcovformula} to calculate the estimated covariance matrix of the residuals($\hat{\Sigma}$)can be one
of the following (see Judge et al., 1955, p.469):
if \texttt{rcovformula}= 0:
\begin{eqnarray*}
\hat{\sigma_{ij}}= \frac{\hat{e_i}\prime\hat{e_j}}{T}
\end{eqnarray*}
if \texttt{rcovformula}= 1 or \texttt{rcovformula}='geomean':
\begin{eqnarray*}
\hat{\sigma_{ij}}= \frac{\hat{e_i}\prime\hat{e_j}}{\sqrt{(T-k_i)\times (T-k_j)}}
\end{eqnarray*}
if \texttt{rcovformula}= 2 or \texttt{rcovformula}='Theil':
\begin{eqnarray*}
\hat{\sigma_{ij}}= \frac{\hat{e_i}\prime\hat{e_j}}{T-k_i-k_j+tr[X_i(X_i\prime X_i)^{-1}X_i\prime X_j(X_j\prime X_j)^{-1}X_j\prime]}
\end{eqnarray*}
if \texttt{rcovformula}= 3 or \texttt{rcovformula}='max':
\begin{eqnarray*}
\hat{\sigma_{ij}}= \frac{\hat{e_i}\prime\hat{e_j}}{T-max(k_i,k_j)}
\end{eqnarray*}
If $i = j$, formula 1, 2, and 3 are equal. All these three formulas yield unbiased estimators
for the diagonal elements of the residual covariance matrix. If $i neq j$, only formula 2
yields an unbiased estimator for the residual covariance matrix, but it is not necessarily
positive semidefinit. Thus, it is doubtful whether formula 2 is really superior to formula 1
(Theil, 1971, p.322).
\subsubsection{Examples}
 Attaching the example dataset:
<<Examples.data>>= 
  data(grunfeld)
@
 Formula:
<<Examples.formula>>=
 formula<-list(mu1=Ige~Fge+Cge,
              mu2=Iw~Fw+Cw)
@ 
Estimating the model using \texttt{sur}:
<<Examples.zelig>>= 
 z.out<-zelig(formula=formula,model="sur",data=grunfeld)
@ 
<<Examples.summary>>=  
summary(z.out)
@ 

Set explanatory variables to their default (mean/mode) values
<<Examples.setx>>= 
 x.out <- setx(z.out)
@ 

Simulate draws from the posterior distribution:
<<Examples.sim>>= 
 s.out <- sim(z.out,x=x.out)
@ 
<<Examples.summary.sim>>=  
summary(s.out)
@ 
\clearpage
\subsubsection{Model}
The basic seemingly unrelated regression model assumes that for each individual
observation $i$ there are $M$ dependent variables ($Y_{ij}, j=1,\ldots,M$) 
each with its own regression equation:
\begin{eqnarray*}
Y_{ij} = X_{ij}'\beta_j + \epsilon_{ij}, \quad  {\rm for} \quad  i=1,\ldots,N \quad {\rm and} \quad  j=1,\ldots,M
\end{eqnarray*}
when $X_{ij}$ is a k-vector of explanatory variables, $\beta_j$
is the coefficients of the explanatory variables,
\begin{itemize}
\item The \emph{stochastic component} is:
\begin{eqnarray*}
\epsilon_{ij}  &  \sim & {\cal N}(0, \sigma_{ij})
\end{eqnarray*}
where within each $j$ equation, $epsilon_{ij}$ is identically
and independently distributed for $i=1,\ldots,M$,
\begin{eqnarray*}
{\rm Var}(\epsilon_{ij})=\sigma_j \quad {\rm and} \quad {\rm Cov}(\epsilon_{ij}, \epsilon_{i\prime j})= 0, \quad {\rm for} \quad i \neq i\prime,\quad {\rm and} \quad j=1,\ldots,M 
\end{eqnarray*}
However, the error terms for the \emph{ith} observation can be correlated across equations
\begin{eqnarray*}
{\rm Cov}(\epsilon_{ij}, \epsilon_{ij\prime})\neq 0, \quad {\rm for} \quad j \neg j\prime, \quad {\rm and} \quad i=1,\ldots,N 
\end{eqnarray*}
\item The \emph{systematic component} is:
\begin{eqnarray*}
\mu_{ij}= E(Y_ij)= X_{ij}\beta_j, \quad {\rm for} \quad  i=1,\ldots,N, \quad {\rm and} \quad j=1,\ldots,M 
\end{eqnarray*}
\end{itemize}
\subsubsection{See Also}
For information about two stage least squares regression, see 
\Sref{2sls} and \texttt{help(2sls)}.
For information about three stage least squares regression, see
\Sref{3sls} and \texttt{help(3sls)}.
\subsubsection{Quantities of Interest}
\subsubsection{Output Values}
The output of each Zelig command contains useful information which you may
view. For example, if you run:
\begin{verbatim}
z.out <- zelig(formula=fml, model = "sur", data)
\end{verbatim}
\noindent then you may examine the available information in \texttt{z.out} by
using \texttt{names(z.out)}, see the draws from the posterior distribution of
the \texttt{coefficients} by using \texttt{z.out\$coefficients}, and view a default
summary of information through \texttt{summary(z.out)}. Other elements
available through the \texttt{\$} operator are listed below:
\begin{itemize}
\item \texttt{rcovest}: residual covariance matrix used for estimation.
\item \texttt{mcelr2}: McElroys R-squared value for the system.
\end{itemize}
\input{systemfit_output}
\begin{itemize}
\item \texttt{maxiter}: maximum number of iterations.
\item \texttt{tol}: tolerance level indicating when to stop the iteration.
\end{itemize}

\subsection* {How to Cite} 

\input{cites/sur}
\input{citeZelig}

\subsection* {See also}
The sur function is adapted from the \texttt{systemfit} library
\citep{HamHen05}.

\bibliographystyle{asa}
\bibliography{gk,gkpubs}
<<afterpkgs, echo=FALSE>>=
 after <- search()
 torm <- setdiff(after,before)
 for (pkg in torm)
 detach(pos=match(pkg,search()))
@
\end{document}

