estimate.auxiliary <- function(x) {
   ar.fit <- ar.ols(x, order.max=3, aic=FALSE, demean=FALSE, intercept=TRUE)
   beta <- c(log(ar.fit$x.intercept), ar.fit$ar, log(sd(ar.fit$resid, na.rm=TRUE)))
   return(beta)
}


lv.indinf <- function(start, data=lynx, s=100) {
    require(numDeriv)
    beta.hat <- estimate.auxiliary(data)
    aux.discrep <- function(theta) {
        aux.from.sims <- replicate(s, estimate.auxiliary(lv.sim(theta.vec=theta)))
        beta.bar <- rowMeans(aux.from.sims)
        sum((beta.hat - beta.bar)^2)
    }
    gradient.discrep <- function(theta) {
        return(grad(func=aux.discrep, x=theta, method="simple"))
    }
   est <- constrOptim(theta=start, f=aux.discrep, grad=gradient.discrep,
                     ui=diag(length(start)), ci=rep(0, length(start)),
                     method="BFGS")
    return(est)
}