R/sigex.mlefit.r

In jlivsey/sigex: Ecce Signum: Multivariate Signal Extraction

#' Fit model to the data using ML estimation
#'
#' @param data.ts A T x N matrix ts object; any missing values
#'			must be encoded with 1i in that entry
#' @param	param  model parameters entered into
#'		a list object with an intuitive structure.
#'    This is an initial specification to
#'			start the nonlinear optimization routines
#' @param	constraint  Matrix of the form [Q , C], with C (constraint.mat)
#'     the matrix of constraints and Q (constraint.vec) the vector
#'     of constraint constants, such that C psi = Q.
#'     Use NULL if there are no constraints
#' @param	mdl  The specified sigex model, a list object
#' @param method  "bfgs" for BFGS, "sann" for simulated annealing,
#'     "cg" for conjugate gradient
#' @param thresh  Pre-parameters theta satisfy |theta|< thresh;
#'     set thresh = Inf if no thresholding is desired
#' @param	hess  A Boolean flag; if true, for BFGS it runs
#'			another round of BFGS to get Hessian matrix
#' @param	whittle  A Boolean flag; if true, uses Whittle likelihood instead of
#'			default Gaussian likelihood
#' @param	debug A Boolean flag; if true, sets the DEBUGGING mode, which prints
#'			psi.last and psi.now to the global field.  Also lik will be printed.
#'			psi.now gives the parameter that crashed the likelihood (if it crashed),
#'			psi.last gives the last good parameter that did not crash the lik.
#' @param maxiter Maximum number of iterations allowed in optimization routines
#'
#' @return 	list with mle and par.est
#'		mle: an object of type outputted by optim
#'		par.est: type param, with the estimated parameters filled in
#' @export
#'

sigex.mlefit <- function(data.ts,param,constraint,mdl,method,thresh=Inf,
                         hess=TRUE,whittle=FALSE,debug=FALSE,maxiter=100)
{

  ##########################################################################
  #
  #	sigex.mlefit
  # 	    Copyright (C) 2017  Tucker McElroy
  #
  #    This program is free software: you can redistribute it and/or modify
  #    it under the terms of the GNU General Public License as published by
  #    the Free Software Foundation, either version 3 of the License, or
  #    (at your option) any later version.
  #
  #    This program is distributed in the hope that it will be useful,
  #    but WITHOUT ANY WARRANTY; without even the implied warranty of
  #    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  #    GNU General Public License for more details.
  #
  #    You should have received a copy of the GNU General Public License
  #    along with this program.  If not, see <https://www.gnu.org/licenses/>.
  #
  ############################################################################

  ################# Documentation #####################################
  #
  #	Purpose: fit model to the data
  #	Background:
  #		param is the name for the model parameters entered into
  #		a list object with a more intuitive structure, whereas
  #		psi refers to a vector of real numbers containing all
  #		hyper-parameters (i.e., reals mapped bijectively to the parameter manifold)
  #	Notes: handles missing values in data.ts, which are indicated by 1i.
  #	Inputs:
  #		data.ts: a T x N matrix ts object; any missing values
  #			must be encoded with 1i in that entry
  #		param: see background; this is an initial specification to
  #			start the nonlinear optimization routines
  #		constraint: matrix of the form [Q , C], with C (constraint.mat)
  #     the matrix of constraints and Q (constraint.vec) the vector
  #     of constraint constants, such that C psi = Q.
  #     Use NULL if there are no constraints
  #		mdl: the specified sigex model, a list object
  #	  method: "bfgs" for BFGS, "sann" for simulated annealing,
  #     "cg" for conjugate gradient
  #   thresh: pre-parameters theta satisfy |theta|< thresh;
  #     set thresh = Inf if no thresholding is desired
  #		hess: a Boolean flag; if true, for BFGS it runs
  #			another round of BFGS to get Hessian matrix
  #		whittle: a Boolean flag; if true, uses Whittle likelihood instead of
  #			default Gaussian likelihood
  #		debug: a Boolean flag; if true, sets the DEBUGGING mode, which prints
  #			psi.last and psi.now to the global field.  Also lik will be printed.
  #			psi.now gives the parameter that crashed the likelihood (if it crashed),
  #			psi.last gives the last good parameter that did not crash the lik.
  #	Outputs:
  #		list with mle and par.est
  #		mle: an object of type outputted by optim
  #		par.est: type param, with the estimated parameters filled in
  #	Requires: sigex.par2psi, sigex.psi2par, sigex.lik, sigex.whittle,
  #   sigex.eta2psi, sigex.psi2eta
  #
  ####################################################################

  fix.lik <- function(eta,constraint,mdl,data.ts,whittle,debug)
  {
    psi <- sigex.eta2psi(eta,constraint)
    # if(debug) {
    #   psi.now <<- psi
    #   message('Current psi:', psi.now)
    # }
    if(whittle)
    {
      out <- sigex.whittle(psi,mdl,data.ts)
    } else
    {
      out <- sigex.lik(psi,mdl,data.ts,debug)
    }
    if(debug)
    {
      psi.last <<- psi
      message('psi = ', paste(psi.last, collapse = ", "))
      message("lik = ", out)
      if(out < lik.best)
      {
        psi.best <<- psi
        lik.best <<- out
      }
    }
    return(out)
  }

  x <- t(data.ts)
  N <- dim(x)[1]
  T <- dim(x)[2]

  par.est <- NULL
  psi <- sigex.par2psi(param,mdl)
  # check: should be zero if constraint holds for initial value
  #	if(length(constraint)>0)
  #	{
  #	  check <- sum((constraint[,-1] %*% psi - constraint[,1])^2)
  #	  print(check)
  #	}
  nueta <- sigex.psi2eta(psi,constraint)
  nu <- nueta[[1]]
  eta <- nueta[[2]]

  # set thresholding to prevent crash (irrelevant for SANN)
  lower.bound <- rep(-thresh,length(eta))
  upper.bound <- rep(thresh,length(eta))

	# initialize best psi
	psi.best <<- sigex.eta2psi(eta,constraint)
	lik.best <<- sigex.lik(psi.best,mdl,data.ts,debug)

	# initial attempt to fit
	if(method=="bfgs") {
	mle <- try(nlminb(eta,fix.lik,constraint=constraint,mdl=mdl,data.ts=data.ts,
	          whittle=whittle,debug=debug,lower=lower.bound,upper=upper.bound,
		        control=list(iter.max=maxiter,eval.max=200)),TRUE)
	}
	if(method=="sann") {
	mle <- optim(eta,fix.lik,constraint=constraint,mdl=mdl,data.ts=data.ts,
	             whittle=whittle,debug=debug,method="SANN",control=list(maxit=maxiter)) }
	if(method=="cg") {
	mle <- optim(eta,fix.lik,constraint=constraint,mdl=mdl,data.ts=data.ts,
	             whittle=whittle,debug=debug,method="CG",
	             control=list(maxit=maxiter,trace=10)) }

  # if the fit was successful, we can report results;
  #	also, we can try to get the hessian - in this case,
  #	we overwrite previous results with L-BFGS-B method
  #	(if not desired, set hess = FALSE)
  if(!inherits(mle, "try-error")) {
    if(hess) {
      eta.est <- mle$par
      mle.hess <- try(optim(eta,fix.lik,constraint=constraint,mdl=mdl,data.ts=data.ts,
                            whittle=whittle,debug=debug,lower=lower.bound,upper=upper.bound,
                            method="L-BFGS-B",hessian=TRUE,control=list(maxit=maxiter)))
      if(!inherits(mle.hess, "try-error")) { mle <- mle.hess }
    }
    eta.est <- mle$par
    psi.est <- sigex.eta2psi(eta.est,constraint)
    par.est <- sigex.psi2par(psi.est,mdl,data.ts) } else mle <- Inf

  return(list(mle,par.est))
}