path.diagram <- function(model, ...){
    UseMethod("path.diagram")
    }

path.diagram.sem <- function(model, out.file, min.rank=NULL, max.rank=NULL, same.rank=NULL,
    variables=model$var.names, parameters=rownames(object$ram), ignore.double=T,
    edge.labels=c("names", "values"), size=c(8,8), node.font=c("Helvetica", 14),
    edge.font=c("Helvetica", 10), rank.direction=c("LR", "TB")){
    if(!missing(out.file)){
        handle <- file(out.file, "w")
        on.exit(close(handle))
        }
        else handle <- stdout()
    edge.labels <- match.arg(edge.labels)
    rank.direction <- match.arg(rank.direction)
    cat(file=handle, paste('digraph "', deparse(substitute(model)),
        '" {\n', sep=""))
    cat(file=handle, paste('  rankdir=', rank.direction, ';\n', sep=""))
    cat(file=handle, paste('  size="',size[1],',',size[2],'";\n', sep=""))
    cat(file=handle, paste('  node [fontname="', node.font[1], 
        '" fontsize=', node.font[2], ' shape=box];\n', sep=""))
    cat(file=handle, paste('  edge [fontname="', edge.font[1],
        '" fontsize=', edge.font[2], '];\n', sep=""))
    cat(file=handle, '  center=1;\n')
    if (!is.null(min.rank)){
        min.rank <- paste('"', min.rank, '"', sep="")
        min.rank <- gsub(',', '" "', gsub(' ', '', min.rank))
        cat(file=handle, paste('  {rank=min ', min.rank, '}\n', sep=""))
        }
    if (!is.null(max.rank)){
        max.rank <- paste('"', max.rank, '"', sep="")
        max.rank <- gsub(',', '" "', gsub(' ', '', max.rank))
        cat(file=handle, paste('  {rank=max ', max.rank,'}\n', sep=""))
        }
    if (!is.null(same.rank)){
        for (s in 1:length(same.rank)){
            same <- paste('"', same.rank[s], '"', sep="")
            same <- gsub(',', '" "', gsub(' ', '', same))
            cat(file=handle, paste('  {rank=same ', same,'}\n', sep=""))
            }
        }
    latent <- variables[-(1:model$n)]
    for (lat in latent){
        cat(file=handle, paste('  "', lat, '" [shape=ellipse]\n', sep=""))
        }
    ram <- model$ram
    ram[names(model$coeff), 5] <- model$coeff
    rownames(ram)[model$fixed] <- ram[model$fixed,5]
    names <- rownames(ram)
    values <- ram[,5]
    heads <- ram[,1]
    to <- ram[,2]
    from <- ram[,3]
    labels <- if (edge.labels == "names") names else values
    direction <- ifelse((heads == 2), ' dir=both', "")
    for (par in 1:nrow(ram)){
        if((!ignore.double) || (heads[par] == 1))
        cat(file=handle, paste('  "', 
            variables[from[par]], '" -> "', variables[to[par]], 
            '" [label="', labels[par], '"', direction[par], '];\n', sep=""))
        }
    cat(file=handle, '}\n')
    }
# last modified 14 April 2001 by J. Fox

print.sem <- function(x, ...) {
    n <- x$n
    t <- x$t
    cat("\n Model Chisquare = ", x$criterion * (x$N - 1), 
        "  Df = ", n*(n + 1)/2 - t, "\n\n")
    print(x$coeff)
    invisible(x)
    }
# last modified 23 Dec 2001 by J. Fox

ram<-function(object, digits=5, startvalues=F){
    old.digits <- options(digits=digits)
    on.exit(options(old.digits))
    var.names <- rownames(object$A)
    ram <- object$ram
    if (!startvalues) colnames(ram) <- c(colnames(ram)[1:4], "estimate")
    par <- object$coeff
    par.names <- rep(" ", nrow(ram))
    t <- object$t
    for (i in 1:t) {
        which.par <- ram[,4] == i
        if (!startvalues)  ram[which.par, 5] <- par[i]
        par.names[which.par] <- names(par)[i]
        }
    par.code <- paste(var.names[ram[,2]], c('<---', '<-->')[ram[,1]],
                    var.names[ram[,3]])
    ram <- data.frame(ram, arrow = par.code)
    print(ram, rowlab=par.names)
    }
# last modified 14 April 2001 by J. Fox

residuals.sem <- function(object, ...){
    object$S - object$C
    }    

standardized.residuals <- function(object, ...){
    UseMethod("standardized.residuals")
    }

standardized.residuals.sem <- function(object, ...){
    res <- residuals(object)
    s <- diag(object$S)
    res/sqrt(outer(s, s))
    }

normalized.residuals <- function(object, ...){
    UseMethod("normalized.residuals")
    }
    
normalized.residuals.sem <- function(object, ...){
    res <- residuals(object)
    N <- object$N
    C <- object$C
    c <- diag(C)
    res/sqrt((outer(c,c) + C^2)/N)
    }
# last modified 23 Dec 2001 by J. Fox

sem <- function(ram, ...){
    if (is.character(ram)) class(ram) <- 'mod'
    UseMethod('sem', ram)
    }

sem.mod <- function (ram, S, N, obs.variables=rownames(S), fixed.x=NULL, debug=F, ...){
    parse.path <- function(path) {
        path.1 <- gsub('-', '', gsub(' ','', path))
        direction <- if (regexpr('<>', path.1) > 0) 2 
            else if (regexpr('<', path.1) > 0) -1
            else if (regexpr('>', path.1) > 0) 1
            else stop(paste('ill-formed path:', path))
        path.1 <- strsplit(path.1, '[<>]')[[1]]
        list(first=path.1[1], second=path.1[length(path.1)], direction=direction)
        }
    if ((!is.matrix(ram)) | ncol(ram) != 3) stop ('ram argument must be a 3-column matrix')
    startvalues <- as.numeric(ram[,3])
    par.names <- ram[,2]
    n.paths <- length(par.names)
    heads <- from <- to <- rep(0, n.paths)
    for (p in 1:n.paths){
        path <- parse.path(ram[p,1])
        heads[p] <- abs(path$direction)
        to[p] <- path$second
        from[p] <- path$first
        if (path$direction == -1) {
            to[p] <- path$first
            from[p] <- path$second
            }
        }
    ram <- matrix(0, p, 5)
    all.vars <- unique(c(to, from))
    latent.vars <- setdiff(all.vars, obs.variables)
    vars <- c(obs.variables, latent.vars)
    pars <- na.omit(unique(par.names))
    ram[,1] <- heads
    ram[,2] <- apply(outer(vars, to, '=='), 2, which)
    ram[,3] <- apply(outer(vars, from, '=='), 2, which)   
    par.nos <- apply(outer(pars, par.names, '=='), 2, which)
    ram[,4] <- unlist(lapply(par.nos, function(x) if (length(x) == 0) 0 else x))
    ram[,5]<- startvalues
    colnames(ram) <- c('heads', 'to', 'from', 'parameter', 'start')
    if (!is.null(fixed.x)) fixed.x <- apply(outer(vars, fixed.x, '=='), 2, which)
    n <- length(obs.variables)
    m <- length(all.vars)
    t <- length(pars)
    if (debug) {
        cat('\n observed variables:\n') 
        print(paste(paste(1:n,':', sep=''), obs.variables, sep=''))
        cat('\n')
        if (m > n){ 
            cat('\n latent variables:\n')
            print(paste(paste((n+1):m,':', sep=''), latent.vars, sep=''))
            cat('\n')
            }
        cat('\n parameters:\n') 
        print(paste(paste(1:t,':', sep=''), pars, sep=''))
        cat('\n\n RAM:\n')
        print(ram)
        }
    sem(ram=ram, S=S, N=N, param.names=pars, var.names=vars, fixed.x=fixed.x)
    }
     

sem.default <- function(ram, S, N, param.names=paste('Param', 1:t, sep=''), 
    var.names=paste('V', 1:m, sep=''), fixed.x=NULL, 
    analytic.gradient=T, heywood=F, warn=F, control=list()){
    is.triangular <- function(X) {
        is.matrix(X) && (nrow(X) == ncol(X)) && 
            (all(0 == X[upper.tri(X)])) || (all(0 == X[lower.tri(X)]))
        }    
    is.symmetric <- function(X) {
        is.matrix(X) && (nrow(X) == ncol(X)) && all(X == t(X))
        }
    if (is.triangular(S)) S <- S + t(S) - diag(diag(S))
    if (!is.symmetric(S)) stop('S must be a square triangular or symmetric matrix')
    con <- list(optim.control=list(),
        optim.method=if (heywood) "L-BFGS-B" else "BFGS", nlm.iterlim=100)
    if ((!is.matrix(ram)) | ncol(ram) != 5 | (!is.numeric(ram)))
        stop ('ram argument must be a 5-column numeric matrix')
    con[names(control)] <- control
    n <- nrow(S)
    observed <- 1:n
    n.fix <- length(fixed.x)
    if (!is.null(fixed.x)){
        for (i in 1:n.fix){
            for (j in 1:i){
                ram <- rbind(ram, c(2, fixed.x[i], fixed.x[j], 
                    0, S[fixed.x[i], fixed.x[j]]))
                }
            }
        }
    m <- max(ram[,2])
    t <- max(ram[,4])
    J <- matrix(0, n, m)
    J[cbind(1:n, observed)]<-1
    par.posn <- unlist(lapply(apply(outer(ram[,4], 1:t, '=='), 2, which), "[", 1))
    colnames(ram)<-c("heads", "to", "from", "parameter", "start value")
    rownames(ram)<-rep("",nrow(ram))
    rownames(ram)[par.posn]<-param.names
    fixed <- ram[,4] == 0
    sel.free <- ram[,4]
    sel.free[fixed] <- 1
    one.head <- ram[,1] == 1
    one.free <- which( (!fixed) & one.head )
    two.free <- which( (!fixed) & (!one.head) )
    arrows.1 <- ram[one.head, c(2,3)]
    arrows.2 <- ram[!one.head, c(2,3)]
    arrows.2t <- ram[!one.head, c(3,2)]
    arrows.1.free <- ram[one.free,c(2,3)]
    arrows.2.free <- ram[two.free,c(2,3)]
    sel.free.1 <- sel.free[one.free]
    sel.free.2 <- sel.free[two.free]
    start <- if (any(is.na(ram[,5][par.posn]))) startvalues(S, ram)
        else ram[,5][par.posn]
    bounds <- rep(-Inf, t)
    bounds[((ram[,1]==2) & (ram[,2]==ram[,3]))[par.posn]] <- if (heywood) .01 else -Inf
    objective.1 <- function(par){
        A <- P <- matrix(0, m, m)
        val <- ifelse (fixed, ram[,5], par[sel.free])
        A[arrows.1] <- val[one.head]
        P[arrows.2t] <- P[arrows.2] <- val[!one.head]
        I.Ainv <- solve(diag(m) - A)
        C <- J %*% I.Ainv %*% P %*% t(I.Ainv) %*% t(J)
        Cinv <- solve(C)
        F <- sum(diag(S %*% Cinv)) + log(det(C))
        F
        }
    objective.2 <- function(par){
        A <- P <- matrix(0, m, m)
        val <- ifelse (fixed, ram[,5], par[sel.free])
        A[arrows.1] <- val[one.head]
        P[arrows.2t] <- P[arrows.2] <- val[!one.head]
        I.Ainv <- solve(diag(m) - A)
        C <- J %*% I.Ainv %*% P %*% t(I.Ainv) %*% t(J)
        Cinv <- solve(C)
        F <- sum(diag(S %*% Cinv)) + log(det(C))
        grad.P <- t(I.Ainv) %*% t(J) %*% Cinv %*% (C - S) %*% Cinv %*% J %*% I.Ainv
        grad.A <- grad.P %*% P %*% t(I.Ainv)
        gradient <- rep(0, m)
        gradient[sel.free.1] <- grad.A[arrows.1.free]
        gradient[sel.free.2] <- grad.P[arrows.2.free]
        attributes(F) <- list(C=C, A=A, P=P, gradient=gradient)
        F
        }
    gradient <- function(par){
        A <- P <- matrix(0, m, m)
        val <- ifelse (fixed, ram[,5], par[sel.free])
        A[arrows.1] <- val[one.head]
        P[arrows.2t] <- P[arrows.2] <- val[!one.head]
        I.Ainv <- solve(diag(m) - A)
        C <- J %*% I.Ainv %*% P %*% t(I.Ainv) %*% t(J)        
        Cinv <- solve(C)
        grad.P <- t(I.Ainv) %*% t(J) %*% Cinv %*% (C - S) %*% Cinv %*% J %*% I.Ainv
        grad.A <- grad.P %*% P %*% t(I.Ainv)
        gradient <- rep(0, m)
        gradient[sel.free.1] <- grad.A[arrows.1.free]
        gradient[sel.free.2] <- grad.P[arrows.2.free]
        gradient
        }
    if (!warn){
        save.warn <- options(warn=-1)
        on.exit(options(save.warn))
        }
    res <- if (con$optim.method == "L-BFGS-B") 
        optim(start, objective.1, method="L-BFGS-B",  
            gr=if (analytic.gradient) gradient,
            control=con$optim.control, lower=bounds)
        else optim(start, objective.1, method=con$optim.method,  
            gr=if (analytic.gradient) gradient,
            control=con$optim.control)
    convergence.1 <- res$convergence
    message.1 <- res$message
    coef.1 <- res$par
    if(res$convergence > 1) res$par <- start
    res <- nlm(if (analytic.gradient) objective.2 else objective.1, 
        res$par, hessian=T, iterlim=con$nlm.iterlim, check=F)
    convergence.2 <- res$code
    if (convergence.2 > 2) {
        coef.2 <- res$estimate
        res$par <- if (convergence.2 > 4) start else coef.2
        res <- nlm(if (!analytic.gradient) objective.2 else objective.1, 
            res$par, hessian=T, iterlim=con$nlm.iterlim, check=F)
        convergence.3 <- res$code
        }
    if (!warn) options(save.warn)
    par <- res$estimate
    names(par) <- param.names
    result <- list()
    result$var.names <- var.names
    obj <- objective.2(par)
    ram[par.posn, 5] <- start
    par.code <- paste(var.names[ram[,2]], c('<---', '<-->')[ram[,1]],
    var.names[ram[,3]])
    result$ram <- ram
    result$coeff <- par
    result$criterion <-  c(obj) - n - log(det(S))
    cov <- (2/(N - 1)) * solve(res$hessian)
    colnames(cov) <- rownames(cov) <- param.names
    result$cov <- cov
    rownames(S) <- colnames(S) <- var.names[observed]
    result$S <- S
    result$J <- J
    C <- attr(obj, "C")
    rownames(C) <- colnames(C) <- var.names[observed]
    result$C <- C
    A <- attr(obj, "A")
    rownames(A) <- colnames(A) <- var.names
    result$A <- A
    P <- attr(obj, "P")
    rownames(P) <- colnames(P) <- var.names
    result$P <- P
    result$n.fix <- n.fix
    result$n <- n
    result$N <- N
    result$m <- m
    result$t <- t
    result$par.posn <- par.posn
    result$convergence.1 <- convergence.1
    result$message.1 <- message.1
    result$coef.1 <- coef.1
    result$convergence.2 <- ultimate <- convergence.2
    if (exists("convergence.3", inherits=F)){
        result$coef.2 <- coef.2
        result$convergence.3 <- ultimate <- convergence.3
        }
    if (convergence.1 > 0) warning('initial optimization DID NOT converge')
    if (convergence.2 > 2) warning('second optimization DID NOT converge')
    if ((convergence.1 > 0) | (convergence.2 > 2)) 
        warning(paste('final optimization', if (ultimate > 2) 'DID NOT' else 'DID','converge'))
    class(result) <- "sem"
    result
    }
# last modified 25 July 2001 by J. Fox

standardized.coefficients <- function(object, digits=5){
    old.digits <- options(digits=digits)
    on.exit(options(old.digits))
    P <- object$P
    A <- object$A
    t <- object$t
    par <- object$coeff
    par.posn <- object$par.posn
    IAinv <- solve(diag(nrow(A)) - A)
    C <- IAinv %*% P %*% t(IAinv)
    ram <- object$ram
    par.names <- rep(' ', nrow(ram))
    for (i in 1:t) {
        which.par <- ram[,4] == i
        ram[which.par, 5] <- par[i]
        par.names[which.par] <- names(par)[i]
        }
    one.head <- ram[,1] == 1
    coeff <- ram[one.head, 5]
    coeff <- coeff *
        sqrt(diag(C[ram[one.head, 3], ram[one.head, 3]])/diag(C[ram[one.head, 2], ram[one.head, 2]]))
    var.names <- rownames(A)   
    par.code <- paste(var.names[ram[one.head,2]], c('<---', '<-->')[ram[one.head,1]],
                    var.names[ram[one.head,3]])
    coeff <- data.frame(par.names[one.head], coeff, par.code)
    names(coeff) <- c(" ", "Std. Estimate", " ")
    print(coeff, rowlab=rep(" ", nrow(coeff)))
    }

std.coef <- function (...){
    standardized.coefficients(...)
    }
# last modified 23 Dec 2001 by J. Fox

startvalues <- function(S, ram){
    n <- nrow(S) 
    observed <- 1:n       
    m <- max(ram[,2])            
    t <- max(ram[,4])   
    s <- sqrt(diag(S))
    R <- S/outer(s,s)
    latent<-(1:m)[-observed]
    par.posn <- unlist(lapply(apply(outer(ram[,4], 1:t, '=='), 2, which), "[", 1))
    one.head <- ram[,1] == 1
    start <- (ram[,5])[par.posn]
    A.pat <-matrix(FALSE, m, m)
    A.pat[ram[one.head, c(2,3)]] <- TRUE
    P.pat <- C <- matrix(0, m, m)
    P.pat[ram[!one.head, c(2,3)]] <- P.pat[ram[!one.head, c(3,2)]] <- 1
    C[observed, observed] <- R
    for (l in latent) {
        indicators <- A.pat[observed, l]
        for (j in observed){
            C[j, l] <- C[l, j] <- if (!any(indicators)) runif(1, .3, .5)
                else {   
                        numerator <- sum(R[j, observed[indicators]])
                        denominator <- sqrt(sum(R[observed[indicators], observed[indicators]]))
                        numerator/denominator
                    }
            }
        }
    for (l in latent){
        for (k in latent){
            C[l, k] <- if (l==k) 1 else {
                                indicators.l <- A.pat[observed, l]
                                indicators.k <- A.pat[observed, k]
                                if ((!any(indicators.l)) | (!any(indicators.k))) runif(1, .3, .5) else {
                                    numerator <- sum(R[observed[indicators.l], observed[indicators.k]])
                                    denominator <- sqrt( sum(R[observed[indicators.l], observed[indicators.l]])
                                        * sum(R[observed[indicators.k], observed[indicators.k]]))
                                    numerator/denominator}
                                    }
            }
        }
    A <- matrix(0, m, m)
    for (j in 1:m){
        ind <- A.pat[j,]
        if (!any(ind)) next
        A[j, ind] <- solve(C[ind, ind]) %*% C[ind, j]
        }
    A[observed,] <- A[observed,]*matrix(s, n, m)
    A[,observed] <- A[,observed]*matrix(s, m, n, byrow=T)
    C[observed,] <- C[observed,]*matrix(s, n, m)
    C[,observed] <- C[,observed]*matrix(s, m, n, byrow=T)
    P <- (diag(m) - A) %*% C %*% t(diag(m) - A)
    P <- P.pat * P
    for (par in 1:t){
        if (!is.na(start[par])) next
        posn <- par.posn[par]
        if (ram[posn, 1] == 1) start[par] <- A[ram[posn, 2], ram[posn, 3]]
            else start[par] <- P[ram[posn, 2], ram[posn, 3]]
        }
    start
    }
# last modified 25 July 2001 by J. Fox

summary.sem <- function(object, digits=5, ...) {
    norm.res <- normalized.residuals(object)
    se <- sqrt(diag(object$cov))
    z <- object$coeff/se
    n.fix <- object$n.fix
    n <- object$n
    t <- object$t
    S <- object$S
    C <- object$C
    N <- object$N
    df <- n*(n + 1)/2 - t - n.fix*(n.fix + 1)/2
    invC <- solve(C)
    CSC <- invC %*% (S - C)
    CSC <- CSC %*% CSC
    CS <- invC %*% S
    CS <- CS %*% CS
    GFI <- 1 - sum(diag(CSC))/sum(diag(CS))
    AGFI <- if (df > 0) 1 - (n*(n + 1)/(2*df))*(1 - GFI)
     else NA
    var.names <- rownames(object$A)
    ram <- object$ram[object$par.posn,]
    par.code <- paste(var.names[ram[,2]], c('<---', '<-->')[ram[,1]],
                    var.names[ram[,3]])
    coeff <- data.frame(object$coeff, se, z, 1 - pnorm(abs(z)), par.code)
    names(coeff) <- c("Estimate", "Std Error", "z value", "Pr(>|z|)", " ")
    row.names(coeff) <- names(object$coeff)
    chisq <- object$criterion * (N - 1)
    BIC <- if (df > 0) chisq - df * log(N*n) else NA
    ans <- list(chisq=chisq, df=df, GFI=GFI, AGFI=AGFI, BIC=BIC, 
        norm.res=norm.res, coeff=coeff, digits=digits)
    class(ans) <- "summary.sem"
    ans
    }
    
print.summary.sem <- function(x, ...){
    old.digits <- options(digits=x$digits)
    on.exit(options(old.digits))
    cat("\n Model Chisquare = ", x$chisq, "  Df = ", x$df, 
        "Pr(>Chisq) =", if (x$df > 0) 1 - pchisq(x$chisq, x$df)
            else NA)
    cat("\n Goodness-of-fit index = ", x$GFI)
    cat("\n Adjusted goodness-of-fit index = ", x$AGFI)
    cat("\n BIC = ", x$BIC, "\n")
    cat("\n Normalized Residuals\n")
    print(summary(as.vector(x$norm.res)))
    cat("\n Parameter Estimates\n")
    print(x$coeff)
    invisible(x)
    }
# Two-Stage Least Squares
#   John Fox

# last modified 21 Nov 2001 by J. Fox

tsls <- function(object, ...){
    UseMethod("tsls")
    }

tsls.default <- function (y, X, Z, names=NULL) {
    n <- length(y)
    p <- ncol(X)
    invZtZ <- solve(crossprod(Z))
    XtZ <- crossprod(X, Z)
    V <- solve(XtZ %*% invZtZ %*% t(XtZ))
    b <- V %*% XtZ %*% invZtZ %*% crossprod(Z, y)
    residuals <- y - X %*% b
    s2 <- sum(residuals^2)/(n - p)
    V <- s2*V
    result<-list()
    result$n <- n
    result$p <- p
    b <- as.vector(b)
    names(b) <- names
    result$coefficients <- b
    rownames(V) <- colnames(V) <- names
    result$V <- V
    result$s <- sqrt(s2)
    result$residuals <- as.vector(residuals)
    result$response <- y
    result$model.matrix <- X
    result$instruments <- Z
    result
    }

    
tsls.formula <- function(model, instruments, data, subset, 
    na.action, contrasts=NULL){
    if (missing(na.action)) 
        na.action <- options()$na.action
    m <- match.call(expand.dots = FALSE)
    if (is.matrix(eval(m$data, sys.frame(sys.parent())))) 
        m$data <- as.data.frame(data)
    response.name <- deparse(model[[2]])
    formula <- as.formula(paste(response.name, '~', 
        deparse(model[[3]]), '+', deparse(instruments[[2]])))
    m$formula <- formula
    m$instruments <- m$model <- m$contrasts <- NULL
    m[[1]] <- as.name("model.frame")
    mf <- eval(m, sys.frame(sys.parent()))
    na.act <- attr(mf, "na.action")
    Z <- model.matrix(instruments, data = mf, contrasts)
    response <- attr(attr(mf, "terms"), "response")
    y <- mf[,response]
    X <- model.matrix(model, data=mf, contrasts)
    result <- tsls(y, X, Z, colnames(X))
    result$response.name <- response.name
    result$formula <- model
    result$instruments <- instruments
    if (!is.null(na.act)) 
        result$na.action <- na.act
    class(result) <- "tsls"
    result
    }

print.tsls <- function(x, ...){
    cat("\nModel Formula: ")
    print(x$formula)
    cat("\nInstruments: ")
    print(x$instruments)
    cat("\nCoefficients:\n")
    print(x$coefficients)
    cat("\n")
    invisible(x)
    }
    
    
summary.tsls <- function(object, digits=4, ...){
    save.digits <- unlist(options(digits=digits))
    on.exit(options(digits=save.digits))
    cat("\n 2SLS Estimates\n")
    cat("\nModel Formula: ")
    print(object$formula)
    cat("\nInstruments: ")
    print(object$instruments)
    cat("\nResiduals:\n")
    print(summary(residuals(object)))
    cat("\n")
    df <- object$n - object$p
    std.errors <- sqrt(diag(object$V))
    b <- object$coefficients
    t <- b/std.errors
    p <- 2*(1 - pt(abs(t), df))
    table <- cbind(b, std.errors, t, p)
    rownames(table) <- names(b)
    colnames(table) <- c("Estimate","Std. Error","t value","Pr(>|t|)")
    print(table)
    cat(paste("\nResidual standard error:", round(object$s, digits),
        "on", df, "degrees of freedom\n\n"))
    }
    
residuals.tsls <- function(object, ...){
    res <- object$residuals
    if (is.null(object$na.action)) 
        res
    else naresid(object$na.action, res)
    }

coefficients.tsls <- function(object, ...){
    object$coefficients
    }
    
fitted.tsls <- function(object, ...){
    yhat <- as.vector(object$model.matrix %*% object$coefficients)
    if (is.null(object$na.action)) 
        yhat
    else napredict(object$na.action, yhat)
    }
    
anova.tsls <- function(object, model.2, s2, dfe, ...){
    if(class(model.2) != "tsls") stop('requires two models of class tsls')
    s2.1 <- object$s^2
    n.1 <- object$n 
    p.1 <- object$p
    dfe.1 <- n.1 - p.1
    s2.2 <- model.2$s^2
    n.2 <- model.2$n
    p.2 <- model.2$p
    dfe.2 <- n.2 - p.2
    SS.1 <- s2.1 * dfe.1
    SS.2 <- s2.2 * dfe.2
    SS <- abs(SS.1 - SS.2)
    Df <- abs(dfe.2 - dfe.1)
    if (missing(s2)){
        s2 <- if (dfe.1 > dfe.2) s2.1 else s2.2
        F <- (SS/Df) / s2
        RSS <- c(SS.1, SS.2)
        Res.Df <- c(dfe.1, dfe.2)
        SS <- c(NA, SS)
        P <- c(NA, 1 - pf(F, Df, min(dfe.1, dfe.2)))
        Df <- c(NA, Df)
        F <- c(NA, F)
        rows <- c("Model 1", "Model 2")
        }
    else{
        F <- (SS/Df) / s2
        RSS <- c(SS.1, SS.2, s2*dfe)
        Res.Df <- c(dfe.1, dfe.2, dfe)
        SS <- c(NA, SS, NA)
        P <- c(NA, 1 - pf(F, Df, min(dfe.1, dfe.2)), NA)
        Df <- c(NA, Df, NA)
        F <- c(NA, F, NA)
        rows <- c("Model 1", "Model 2", "Error")
        }
    table <- data.frame(Res.Df, RSS, Df, SS, F, P)
    head.1 <- paste("Model 1: ",format(object$formula), "  Instruments:", 
        format(object$instruments))
    head.2 <- paste("Model 2: ",format(model.2$formula), "  Instruments:", 
        format(model.2$instruments))
    names(table) <- c("Res.Df", "RSS", "Df", "Sum of Sq", "F", "Pr(>F)")
    row.names(table) <- rows
    structure(table, heading = c("Analysis of Variance", "", head.1, head.2, ""), 
        class = c("anova", "data.frame"))
    }
